feature extraction for a pretrained model in keras
The features variable contains the outputs of the final convolutional layers of your network. The final convolutional layer of VGG16 outputs 512 7x7 feature maps. All you need to do in order to use these features in a logistic regression model (or any other model) is reshape it to a 2D tensor, as you say. 

reshaped_features = features.reshape(100, 512*7*7)


This will flatten out the feature maps to a long one-dimensional vector for each instance.

Tip: If you can't be bothered working out the actual dimensions for your reshape, you can replace 512*7*7 with -1 and numpy will figure out how big the final dimension should be.