texts: Selecting certain road length around points in QGIS

texts:
I have a dataset with road segments across a country, and their length in km.
I also have the boundaries of each province within the country.
I have calculated the center point of each province.
What I would like to do now in QGIS is to select a contiguous % of road lengths per province, and have the selection radiate out from the center point of each province.
Here's a screenshot of what I'd like to achieve: rather than randomly selecting 10% of the roads (shown in green and orange) I have the red dot, and want to just highlight 10% of roads (by length) around the dot (as shown in blue).


texts: Is GIS industry always ruled by two options only, namely Esri products and OSGeo projects?

texts:
I always wonder, why in GIS industry there are merely two major options available. Although there are many other proprietary software products for GIS, their market share is well below ESRI products (here, here and here although quite back dated). My opinion became stronger after joining GIS SE. Most of the questions, as well as the leading users (with highest reputation), are either for/interested in ESRI line of products or FOSS (of course apart from some questions related to data and theory).

However this not the case in other streams of information technology, databases and programming languages in broad sense. I personally feel one reason could be the lack of interest by the software giants like Microsoft, Google or Oracle. While the first two have WebMaps, these are not WebGIS, and the spatial extension of Oracle is not fully GIS too.

What are the major reasons for such a trend? There is nice article somewhat related to this, " The Five Major Roadblocks to GIS Gaining Corporate Market Share". But this is not the answer I am looking for.

Most of the answers seem to be trying to prove open source tools are better than Esri products! But that was not the question. Because it is difficult to compare a matured but expensive tool with a bunch of comparatively newer open source tools with great usability. My question is why no industry standard product from any established software giants?


texts: ArcGIS dashboard Arcade Using multiple If/else if statements

texts:
I try to write an Arcade expression in dashboard indicator, the expression has 4 if statements but only the first two statements are works the code is below:
############################################

var`enter code here`ar kmillion = 10000000;
var million = 1000000;
var k100 = 100000;

var name_km = Left($datapoint["sum_VALUE"],2)
var name_m  = Left($datapoint["sum_VALUE"],1)
var name_k100 = Left($datapoint["sum_VALUE"],3)
var name_k99 =  Left($datapoint["sum_VALUE"],2)

if ($datapoint["sum_VALUE"] >= kmillion) {
    return {
        topText: 'اجمالي عدد التلاميذ',
        topTextColor: '#a80000',
        topTextOutlineColor: '',
        topTextMaxSize: 'small',
        middleText: name_km,
        middleTextColor: '',
        middleTextOutlineColor: '',
        middleTextMaxSize: 'large',
        bottomText: 'مليون تلميذ',
        iconName:'icon1',
        iconAlign:'left',
        iconColor:'',
        iconOutlineColor:''

    }
}

else if ($datapoint["sum_VALUE"] >= million < kmillion  ) {
    return{
        topText: 'اجمالي عدد التلاميذ',
        topTextColor: '#a80000',
        topTextOutlineColor: '',
        topTextMaxSize: 'small',
        middleText: name_m,
        middleTextColor: '',
        middleTextOutlineColor: '',
        middleTextMaxSize: 'large',
        bottomText: '2مليون تلميذ',
        iconName:'icon1',
        iconAlign:'left',
        iconColor:'',
        iconOutlineColor:''

    }
}

else {
    return {
        topText: 'اجمالي عدد التلاميذ',
        topTextColor: '#a80000',
        topTextOutlineColor: '',
        topTextMaxSize: 'small',
        middleText: name_k100,
        middleTextColor: '',
        middleTextOutlineColor: '',
        middleTextMaxSize: 'large',
        bottomText: 'ألف تلميذ',
        iconName:'icon1',
        iconAlign:'left',
        iconColor:'',
        iconOutlineColor:''

    }
}


texts: OpenStreetMap and Python recursively find all the bus routes in a city

texts:
I want to retrieve all the bus routes (as a collection of lines) from openstreetmap for a given city.

I'm currently using the overpy module (v0.4) to access the Overpass API and this query:

QUERY = """
[out:json];
area[name="%(city)s"]->.a;
(
  node["route"="%(type)s"](area.a);
  way["route"="%(type)s"](area.a);
  relation["route"="%(type)s"](area.a);
);
out body;
>;
out skel;
"""

where city is the city i'm interested in, and type could either be bus or tram.

The problem is, this does not work often :( If I'm looking for the bus routes for a very well known city like Chicago it returns an empty set, while other times the result set is partial (like for Philadelphia for example); very rarely it is able to correctly select the city and grab all the routes.

If I use the same query (without the area selector) in http://overpass-turbo.eu/ when focused on Chicago it is able to return the routes as expected; so I suspect there's something wrong with the way I'm using area

Do you have any idea how can I achieve this?


texts: I have created slope and aspect from a DEM layer, I reclassed both slope and aspect into categories (for example slope is now a value of 1-6) and turned them into polygon layers where each polygon has a different slope value etc. At the moment they are separate so I have polygons for aspect and polygons for slopes and I want to end up with one polygon layer where each polygon has a slope value and an aspect value, does anyone know how I might do this?


texts:
I have QGIS V 2.18.18

I am practicing making 3D maps and decided to do a simple population density 3D map for the states and territories of India.
I have a shapefile that contains 36 features(States and territories) with their respective densities.

I used the 2.5D option in the layer styling, and chose pop_density/1000 as the height.

Now it is rendering, but it is taking a really really long time to render the map.

It has been 15 minutes and less than 1/3rd of the map has been rendered.

I have an Intel Corei7 2.9 GhZ processor, 32 gb of memory, and and NVIDIA quadro 1200M on my laptop. So I do not think there could be problem with my computer configuration.

Should I be checking some configuration, or changing any settings for it to render faster?


texts: Is it possible to figure out which "Path"(PointCollection) is clicked on within a PolyLine?

texts:
I am using ArcGIS WPF's Polyline to plot several lines (PointCollections). Adding a Path(PointCollection) looks something like this.

PointCollection collection = new PointCollection();
PolyLine myLine = new PolyLine();

collection.Add(firstMapPoint);
collection.Add(secondMapPoint);

myLine.Paths.Add(collection);

PolyLine is a Geometry in which I can set to a Graphic and plot on a map. Once that Graphic is clicked it of course is clicking the entire graphic. I am visually showing several Paths(PointCollections). Is it possible to figure out which Path is being clicked on from the user?

Note

I have a solution now where I just have 1 PointCollection per PolyLine. So I am adding one PolyLine/Graphic for each line. I would rather not do this if I can avoid it.


texts: Georeferencing aerial photos when only centroid is known using ArcGIS for Desktop or ERDAS Imagine?

texts:
I have heaps of aerial photos that need to be georeferenced.

I have access to ArcMap and ERDAS.

All I am given is an excel sheet with the centroid coordinates of each image.

Each photo covers an 18m x 24m region.

The photos are in a very remote region with no roads or other structures to georeference them to.

Is it possible to do this?


texts: Adding geometry to QgsVectorLayer with lon/lat coordinates

texts:
I'm creating a custom Qt app that uses pyqgis for showing data on maps. I have successfully created a map canvas that shows a world map from a WMS server.

Now I'd like to overlay data on this map. However, it seems that my vector layer interprets all coordinates in meters instead of lon/lat.

Here's more or less what I do to setup my data layer:

my_layer = QgsVectorLayer('Point', 'My Layer', 'memory')
my_layer.setCrs(QgsCoordinateReferenceSystem(4326))
pr = my_layer.dataProvider()
pr.addAttributes([QgsField('name', QVariant.String),
                  QgsField('value', QVariant.Double)])

# add a point at lon=7.8, lat=47.5
feature = QgsFeature()
feature.setGeometry(QgsGeometry.fromPoint(QgsPoint(7.8, 47.5)))
feature.setAttributes(['some name', 0.0])
pr.addFeatures([feature])

# ... add layer to registry and canvas

The point shows up nicely but in the wrong location. How can I get qgis to interpret my coordinates in lon/lat?


texts: Creating new raster with minimum (or maximum) value betwen two rasters in ArcMap?

texts:
I am using ArcMap 10.3.

I have two temperature rasters from two months. I need to create a new raster from these, but the new raster should have in each cell the maximun value betwen the initial raster.
Example.

Raster 1     Raster 2       Output I need
15             18                  18
19             11                  19


texts: Geotools Label location issue with custom Symbol

texts:
We have a Feature Layer having styling with a symbol using some font create by us.
But we need to apply Labels on this following is the result we expected from this:-

and following is the result we have achieve with the AnchorPoint ,Displacement and Rotation :-

Code For Label Style:-

Rule rule = styleFactory.createRule();
TextSymbolizer textsym = Styling.styleFactory.createTextSymbolizer();
textsym.setLabel(Styling.filterFactory.property("ColumnName"));
StyleBuilder sb = new StyleBuilder();
textsym.setFill(sb.createFill(Color.black));
AnchorPoint anchorPoint = sb.createAnchorPoint(sb.attributeExpression("the_geom"),
        sb.attributeExpression("the_geom"));
Displacement displacement = sb.createDisplacement(filterFactory.literal(10), filterFactory.literal(4));
PointPlacement pointPlacement = sb.createPointPlacement(anchorPoint,displacement,filterFactory.property("RotationColumn"));
textsym.setLabelPlacement(pointPlacement);
rule.symbolizers().add(textsym);


texts: QGIS 2.14.14 Essen LRS Error "TypeError: object of type 'NoneType' has no len()"

texts:
I have recently installed the Essen version of QGIS and I am having difficulty running a query I have ran without problem on previous versions of QGIS.

I want to provide 'measurement' data to a polyline using a points layer by calibrating with LRS.

I receive the following Python Error

  TypeError: object of type 'NoneType' has no len()

Is this a bug in the latest version or is there something simple I am overlooking?


texts:
In the view menu, there's an option called "Toggle Map Only." I selected that option to see what it looks like. This is what the entire QGIS window looks like in Map Only view:

As you can see, this option turns off all of the menus, panels and other visible controls. The only way to disable it is to use the hotkey combo. Of course, if you didn't memorize the hotkey combo, you can't find it again, because you don't have access to the View menu anymore.

How to un-toggle "Map Only" view mode?


texts: Using rasterio to crop image using pixel coordinates instead of geographic coordinates

texts:
I'm trying to extract a series of random patches of image from a larger satellite image (from worldview3). I want to extract patches of uniform size (say 512x512 for example). If it were a normal image, any number of libraries could do this easily. But I want to use rasterio to retain the geographic information of the image patch.
So far the code I've written identifies the upper left pixel coordinate using
image.bounds

and using the resulting bounding box to get the height and the width. Then I use numpy to randomly generate a new random coordinate within the image extent:
new_col = np.random.randint(min(bounds.left, bounds.right), max(bounds.left, bounds.right)-(size+1))
new_row = np.random.randint(min(bounds.top, bounds.bottom), max(bounds.top, bounds.bottom)-(size+1))

And then I subtract the patch size (512) from the new row and column and use those new coordinates as the minx, miny, maxx, maxy and crop from there using
mask(image, shapes=coords, crop=True)

With a non-georeferenced image where the upper left is (0,0) and the lower right is (M,N) this works flawlessly. Similarly for a NAIP image this seems to work. But with the worldview3 image the size is not uniform. I'll get an image size like 700x2000 for example, where I want it to be 512x512. My thought was that the NAIP image I tested with was in UTM, but that the worldview3 image was in lat/long so that subtracting the patch size from the random row and column didn't translate to a uniform pixel size. I still think the error has to do with lat/long but after reprojecting the image to UTM the problem persists.
So, is there any way I can crop the worldview3 image using rasterio but using pixel coordinates instead of geographic coordinates so that I can crop a uniform image size from the larger image but still retain the geographic information of the cropped image patch?


texts: Distance from points to borders of polygons in QGIS big data set

texts:
In QGIS I have two shapefiles: one point layer representing power plants (654 elements) and one polygon layer representing cities and villages connected to gas grid (1348 elements).

For every power plant I need to find all the areas, which are maximum in 10 km distance from the plant. What I can imagine in QGIS I need to find all the areas to whose border the distance from the point is smaller or equal 10 km.

As a result I want to obtain a table with the list of power plants, selected areas and their distances.

How can I make it in QGIS? I am familiar with Python, but never worked with PostGIS or GRASS.


texts: Problems when converting KML files to GPX

texts:
Looking at the KML and the GPX file that I created (in QGIS through 'Save as'), in QGIS it looks as if the lines are all right (picture 1), but when I open it in another type of GPS program like bikegpx.com, the lines are all weird..

Can anyone tell me, how to fix this problem?


texts: My KMZ works for older versions of GE, but not the newer ones

texts:
I have a KMZ file which is a time animation. It used to work in the older GE, such as version 5 (I think). I can load the KMZ file in the newer GE, but do not see the animation control bar to play it. Something needs to be changed in the KMZ file, but cannot figure out. Here is the beginning and end of the KML file:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.0">
<Folder>
  <name maxLines="1">ADCIRC_pbl</name>
  <Snippet maxLines="1"><![CDATA[]]>
  </Snippet>
  <Folder><name><![CDATA[ahframe_legend]]></name><visibility>1</visibility><ScreenOverlay><visibility>1</visibility><overlayXY x="1" y="1" xunits="fraction" yunits="fraction" /><screenXY x="1" y="1" xunits="fraction" yunits="fraction" /><Icon><href>ahframe_legend.png</href></Icon></ScreenOverlay></Folder>
<Folder>
<name maxLines="1">ahframe_1: Image</name>
  <Snippet maxLines="1"><![CDATA[]]>
  </Snippet>
  <GroundOverlay>
    <TimeSpan>
        <begin>2005-08-28T23:00:00z</begin>
        <end>2005-08-28T23:30:00z</end>
    </TimeSpan>
    <name maxLines="1">ahframe_1</name>
    <Icon>
      <href>ahframe_1.png</href>
      <viewBoundScale>1.0</viewBoundScale>
    </Icon>
    <LatLonBox>
      <north>30.46</north>
      <south>29.44</south>
      <east>-88.3531</east>
      <west>-89.6952</west>
       <rotation>0</rotation>
     </LatLonBox>
   </GroundOverlay>
 </Folder>
<Folder>
<name maxLines="1">ahframe_2: Image</name>
  <Snippet maxLines="1"><![CDATA[]]>
  </Snippet>
  <GroundOverlay>
    <TimeSpan>
        <begin>2005-08-28T23:30:00z</begin>
        <end>2005-08-29T00:00:00z</end>
    </TimeSpan>
    <name maxLines="1">ahframe_2</name>
    <Icon>
      <href>ahframe_2.png</href>
      <viewBoundScale>1.0</viewBoundScale>
    </Icon>
    <LatLonBox>
      <north>30.46</north>
      <south>29.44</south>
      <east>-88.3531</east>
      <west>-89.6952</west>
       <rotation>0</rotation>
     </LatLonBox>
   </GroundOverlay>
 </Folder>
<Folder>

The KML file ends like this:

<name maxLines="1">ahframe_77: Image</name>
  <Snippet maxLines="1"><![CDATA[]]>
  </Snippet>
  <GroundOverlay>
    <TimeSpan>
        <begin>2005-08-30T13:00:00z</begin>
        <end>2005-08-30T13:30:00z</end>
    </TimeSpan>
    <name maxLines="1">ahframe_77</name>
    <Icon>
      <href>ahframe_77.png</href>
      <viewBoundScale>1.0</viewBoundScale>
    </Icon>
    <LatLonBox>
      <north>30.46</north>
      <south>29.44</south>
      <east>-88.3531</east>
      <west>-89.6952</west>
       <rotation>0</rotation>
     </LatLonBox>
   </GroundOverlay>
 </Folder>
 </Folder>
 </kml>


texts: ArcGIS Pro Map Series - prevent overlapping inset map

texts: I'm using Map Series in ArcGIS Pro to produce layouts for around 1,000 small area boundaries. The layout includes an inset map in the top right corner showing the location of the small area in the main map at a larger scale. However, in some of the layouts, the inset map overlaps with a part of the small area. Is there an automated method that places the inset map in each layout in the Map Series in an appropriate corner (top right, top left, bottom right or bottom left) to ensure that it doesn't overlap with the small area boundary in the main map?


texts: How to handle A column was specified that does not exist exception. and still put a value where the column should be

texts:
ArcGIS 10.4

I have a python (arcpy) script that is pulling values from a input list of fields for every featureclass in a geaodatabase (gdb) using a SearchCursor and is putting those values in a csv.  The script works perfect for all featureclasses that have all the fields in the list.  The script below is a thought process of how to capture all the field names in a list and do something if they are all present.  This works but hits a stop if a field does not exist.

I've handled the " A column specified does not exist" using try/except exception in a different script that does not use the listFields so the script continues to run, but if even a single field does not exist in that featureclass the entire featureclass is skipped and is not present in the csv.

I need all records from all featureclasses in the csv even if a field does not exist in that featureclass.  So something like if fieldname[0] not in fieldList put "field does not exist" where otherwise that fields value would be.  I need this for every field in the list. The csv will have featureclassname, field1, field2, field3..etc.

import arcpy
import csv
import operator

arcpy.env.workspace = r"path\to\my\gdb"
newCSV = r"path\to\new\csv"
valueList = []

# Use the ListDatasets function to return a list of all feature datasets in the gdb:
datasetList = arcpy.ListDatasets("*", "Feature")

with open(newCSV,'wb') as csvfile:
    for dataset in datasetList:

# Use the ListFeatureClasses function to return a list of all fc's in the gdb:
    fcList = arcpy.ListFeatureClasses("*","",dataset)
        for fc in fcList:

        fieldname = ['field1', 'field2', 'field3']
        fieldList = [f.name for f in arcpy.ListFields(fc)]

        if fieldname[0] in fieldList and fieldname[1] in fieldList and fieldname[2] in fieldList:

            with arcpy.da.SearchCursor(fc, fieldname) as cursor:
                row = cursor.next()

                while row:
                    valueList.append([[fc,row[0],row[1],row[2]]])

                    # #merge lists of lists
                    flattenList =  reduce(operator.concat,valueList)

                    writer = csv.writer(csvfile)
                    writer.writerows(flattenList)

                    row = cursor.next()
                    valueList = []


texts: Conventions for returning GeoJSON via API

texts:
I have a queryable API that returns one or more GeoJSON FeatureCollections. Currently the response looks something like this:

{
  "success": {
    "data": {
      "layer1": {"type": "FeatureCollection"....},
      "layer2": {"type": "FeatureCollection"....}
    }
}

However, while trying to view the results in QGIS I realized this might not be the best approach (as QGIS doesn't allow you to specify an access key when adding a GeoJSON layer via a URI). I think this is generally a sensible design, except that I can't find any other examples of returning GeoJSON as part of a larger JSON structure.

Is there anything wrong with a GeoJSON feature (specifically a FeatureCollection, but really any type) being part of another "document"?

Though GeoJSON is JSON, I can also see a strong argument that if a user requests a resource as GeoJSON, it is rational to expect a valid GeoJSON document in return.


texts: Layers of Shapefiles are not overlapping

texts:
I have analysed Spacial Poisson of the data from Nepal by using SaTscan and tried to join all the coordinates BY QGIS but it shows map and the points are in different locations.

I have already checked the X Y code is OK and matched with Nepal Data.


texts: How to make symbols visible in legend in Layout?

texts:
I have followed the steps for creating proportional symbol and multivariate analysis (QGIS docs).
It looks like expected in the canvas and in the Legend Items-window in "Item properties"-tab in the Layout.
However, the symbols are not visible in the legend in the map.
How may I make the symbols visible also here?
I'm using QGIS ver. 3.16.0.


texts: First click layer selector doesn't work in CartoDB?

texts:
I'm working on an interactive map using CartoDB built of Leaflet.js. I added a layer selector which works fine after the first click. The first click just doesn't work, if I click the selected layer (Layer 2) it reloads it, but click any other layer and the data is removed from the map, until a second click. Does anyone know the solution for the first click not responding?

The html of the layer selector:

LAYER SELECTOR

Layer 1

Layer 2

Route 22

And the js detecting the click and creating the selector:

var vizjson = 'https://thomvandegeer.cartodb.com/api/v2/viz/a8fc004a-b86b-11e5-b7ca-0ecd1babdde5/viz.json';
cartodb.createLayer(map_object, vizjson)
.addTo(map_object)
.done(function(layer) {
    $(".mt-layer-select").on('click', function(e) {
        var num = +$(e.target).attr('data');
        createSelector(layer,num,$(e.target));
    })
    $(".route-22").on('click', function(){
        map_object.setView([52.3665, 4.897], 17, {animation: true});
    })
    $(".reset-view").on('click', function(){
        map_object.setView([52.37, 4.9], 15, {animation: true});
    })
});

function createSelector(layer,num) {
    for (var i = 0; i < layer.getSubLayerCount(); i++) {
        if (i === num) {
            layer.getSubLayer(i).show();
        } else {
            layer.getSubLayer(i).hide();
        }
    }
}

And it's online here.


texts: How to calculate values through nested grids in QGIS/GeoPandas

texts:
Problem: I have a series of nested grids. Example is below. I want to perform calculations for the smallest grid (grid with the smallest cells, green here) based on the values in the upper cells without performing an intersection method due to performance issues (my final grid has 4^10 cells).
Example:
Assume, that we have 4 grids. The biggest one, red, has only 1 cell with the value 'A' in it. At the next level we have blue grid with 4 cells and values 'X, Y, B, Z' in it, etc. At the lowest level I have pink grid with 64 cells.
I want to calculate the list of upper values for the cells at the smallest grid. So, for the selected pink cell it will be: [A, B, C, D] because [A is the the value in the upper cell at the 1st grid, B is the value in the upper call at the 2nd grid and so on]

Solution ideas: As I've said, the idea here is not to perform intersection method due to its obvious exponential performance issues for the grid with a HUGE number of cells (4x2^10 at 10th level in my case, so it has to intersect each cell with 4x2^9 cells at the 9th level, 4x2^8 cells at the 8th, ...)
My idea was to somehow implement cell IDs based on their place in the grid. Walking from upper left to lower right corner and setting 1,2,3 or 4 value I get IDs making possible to split it based on the level I want to get value from. The problem here is to inherit values from the upper cells.

I can use GeoPandas, QGIS or PyQGIS if you write down the tutorial because I am not very familiar with PyQGIS syntax.


texts: How do I check which version of GDAL I have and how many?

texts:
I have installed in my computer the following software:

PosgresSQL 9.3.2.3
PostGIS 2.1.1.1
QGIS 2.0.1 Dufour (this one installed through the OSGeo4w package)
QGIS 1.7.0 Wroclaw

As far as I know all of them somehow operate, or use, or access the GDAL library.

Are they independent versions of GDAL therefore being used depending on which one of the above software I am using? For instance, If I open the OSGeo4W Shell, just after loading the batch it displays:

GDAL 1.10.1, released 2015/02/10

So, how can I check the other GDAL versions and make comparisons?

Is there any way of centralize all of this stuff and force the apps to use the same centralized folder?


texts: I have created a model that allows the user to delete and edit certain fields in a feature class (Simple model for first time users of ArcGIS so they don't mess up data). The first step of the model is to copy the feature class so we have an original copy to the data. I am curious if I can edit the output name of the copy feature class tool. Currently the output is : test_CopyFeatures. I would like Test_Original.


texts: Unable to edit the "id" field in in a qgis shapefile layer (but everything else seems to work)

texts:
I am trying to create a polygon shapefile layer in QGIS (2.12.1).  Because I will be creating many polygons I'm not going to be able to maintain a unique id number manually (without keeping track of it on a piece of paper or something, which is exactly the sort of thing these fancy computers were supposed to replace).

I'm trying to follow these instructions to automatically assign an id after the fact, but every attempt I make to edit the "id" column fails to persist.  As soon as I save the edits (or just try sorting the attribute table on the id column or something) everything reverts to NULL.

I see the same behavior in 2.12.0 as well.

Is this supposed to work?  Am I missing something obvious?


texts: Only Capability appearing is Feature Access

texts: I was wanting to know if anyone could explain why the only capability showing when I try to publish an ArcGIS Service is Feature Access (Picture Below). I am needing to create a service that will allow the user to update, add, and delete attributes from a feature layer. I have the database registered, but the only Feature Access capability I am able to access is the ability to query. Any help would be appreciated.


texts: How to get full 16 bits range using RasterSymbolizer and ColorMap ramp?

texts:
I am attempting to get the full 16 bits values of my DEM (stored as TIFF files) rendered in WMS tiles (png format).

For this I use GeoServer 2.5 and a RasterSymbolizer with a ColorMap of type "ramp".
Then I used the color "#000000" for the value 0 and the color "#003E80" for 15000. My understanding was that GeoServer will interpolate as such:
1 -> #000001
2 -> #000002
256 -> #000100
15000 -> #003E80
etc...
Instead it looks like it is taking the 16 bits values and divide it by 255 and apply the coefficient to each color. This in effect reduces my height resolution by a big amount!

Is there a way to have the behavior I expected using a ColorMap (or any other sld operator)?

Thanks.


texts: QGIS & WMS: Different aspect ratio between width/height and bbox

texts:
I am experiencing a small screen space offset when using a WMS layer (origin is GeoServer 2.1.4). Looking at the requests from QGIS (1.7.4) to my GeoServer I realised that the aspect ratio between the WIDTH/HEIGHT and the BBOX is different. Calculating the aspect difference back into screen space gives a difference of ~1 pixel - I wonder if it could be the reason for my offset?

Due to devices with non-square pixels it makes sense that the aspect ratio can differ.

But what is the reason for QGIS not enforcing equality between the aspect ratios?


texts: How to move a Point perpendicular from nearest Line?

texts:
I want to move Point perpendicular from nearest line or polygon, and the question is how to do that?

example, I have points (house number) near roadbed, and I don't want them to intersect roadbed (buffered street line)

I have used near analysis to get street line angle to the points for correct labeling, so the angle of the labels is really good but they are to close the road.

So I want to move the points perpendicular away from street line. is there any way?


texts: Change color of top of interstate shield in ArcGIS Pro

texts:
I would like to know how to change the color of the top of an interstate shield text symbol in ArcGIS Pro. The default interstate shield symbol is a blue shield with a red top. In the text symbol properties, it is easy to change the bottom of the shield from blue to another color (see screenshot). How can I change the top of the shield? I am able to do this in ArcMap but can't identify how to do it in ArcGIS Pro (2.1.0).

I want interstates labeled, but want the labels to be more muted, so I am trying to make them gray instead of colored.


texts: I installed version 3.4.0 and repeatedly searched for the QTILE3(QTiles and QTile , i try both of them) plugin, but failed. Why do not I find this plugin?can anyone help me?


texts: Select only records with region objects in MapBasic

texts:
I am writing a program that finds overlapping objects. I can't perform the AreaOverlap() function on objects that don't enclose an area (e.g. text, points) as the result would be meaningless.

Is there a way to select only records that have objects that are regions?

I have a function that checks whether an object is a region (see below) so I could loop through the objects and check them one-by-one but it would be more efficient to do this as part of a select statement.

My example function:

Function IsObjectRegion(ObjAlias As Object) As Logical
  Do Case ObjectInfo(ObjAlias, OBJ_INFO_TYPE)
    Case OBJ_TYPE_REGION,
         OBJ_TYPE_RECT,
         OBJ_TYPE_ROUNDRECT,
         OBJ_TYPE_ELLIPSE
      IsObjectRegion = True
    Case Else
      IsObjectRegion = False
  End Case
End Function

Update: my final, working code using Peter's method and my function:

Include "mapbasic.def"

Declare Sub main
Declare Function IsObjectRegion(ByVal ObjAlias As Object) As Logical

Function IsObjectRegion(ByVal ObjAlias As Object) As Logical
  Do Case ObjectInfo(ObjAlias, OBJ_INFO_TYPE)
    Case OBJ_TYPE_REGION,
         OBJ_TYPE_RECT,
         OBJ_TYPE_ROUNDRECT,
         OBJ_TYPE_ELLIPSE
      IsObjectRegion = True
    Case Else
      IsObjectRegion = False
  End Case
End Function

Sub Main
  Select * from TargetTable Where obj into temp1 noselect
  Select * from temp where IsObjectRegion(OBJ) into results noselect
  Browse * from results
  Close table temp
End Sub


texts:
I am doing a research project for which I need to find an effective way for drawing lines between a number of points along a square. I need there to be a line between every point on the square in order to simulate the ability to cross the square whichever way a pedestrian pleases.

Ideally the process should generate a new line layer.

The end result needs to look something like this:

I have taken a look at previous posts on the subject, however the plugins that they recommend - PointsToPaths, Points2One - don't seem to be available in QGIS anymore, and the Points2One one doesn't work.


texts: GRASS v.to.rast outputs a bounding box of the geographic region

texts:
I'm trying to rasterize a vector layer in GRASS after setting the geographic region to the vector layer extent. I also used v.category to add categories to my vector layer. The command I'm using to rasterize is:
v.to.rast -d --overwrite input=nyc_census_blocks@may type=boundary,centroid,area output=raster_map@may use=cat

However, the output I'm getting is just a yellow rectangle of the geographic region I set, not the shape of my vector layer.
The region is:
g.region -p
projection: 3 (Latitude-Longitude)
zone:       0
datum:      wgs84
ellipsoid:  wgs84
north:      40:54:55.917994N
south:      40:29:46.015423N
west:       74:15:20.12891W
east:       73:42:00.03263W
nsres:      0:25:09.902571
ewres:      0:33:20.09628
rows:       1
cols:       1
cells:      1

Description of vector layer:
v.info nyc_census_blocks

+----------------------------------------------------------------------------+
 | Name:            nyc_census_blocks                                         |
 | Mapset:          may                                                       |
 | Location:        newLocation                                               |
 | Database:        /Users/may                                                |
 | Title:                                                                     |
 | Map scale:       1:1                                                       |
 | Name of creator: may                                                       |
 | Organization:                                                              |
 | Source date:     Thu Mar  4 11:32:24 2021                                  |
 | Timestamp (first layer): none                                              |
 |----------------------------------------------------------------------------|
 | Map format:      native                                                    |
 |----------------------------------------------------------------------------|
 |   Type of map: vector (level: 2)                                           |
 |                                                                            |
 |   Number of points:       0               Number of centroids:  39250      |
 |   Number of lines:        0               Number of boundaries: 91570      |
 |   Number of areas:        39250           Number of islands:    107        |
 |                                                                            |
 |   Map is 3D:              No                                               |
 |   Number of dblinks:      1                                                |
 |                                                                            |
 |   Projection: Latitude-Longitude                                           |
 |                                                                            |
 |               N:  40:54:55.917994N    S:  40:29:46.015423N                 |
 |               E:   73:42:00.03263W    W:   74:15:20.12891W                 |
 |                                                                            |
 |   Digitization threshold: 0                                                |
 |   Comment:                                                                 |
 |                                                                            |
 +----------------------------------------------------------------------------+


texts: QGIS Spatial Index uses bounding boxes; how to get actual nearest line to a known point

texts:
I have looked at lots of answers referring to this same issue, and tried using QGIS Spatial Indexes and iterating through the line features in the spatial index to calculate the distance between a point and a line, but pt.distance requires either pt.distance(float, float) or pt.distance(QgsPoint). pt.distance won't accept nnfeature.geometry() if that geometry is a line. Is there a way to refer to line as points that are 1m apart?

# Find a nearest index entry to myPoint (QgsPoint)
nearestindexid = spacialIndex.nearestNeighbor(pt, 1)[0]
# Get the feature in the layer (self.myvectorlayer) that is
# being indexed by self.spIndex
nnfeature = featureID_Layer.getFeatures(QgsFeatureRequest(nearestindexid)).next()
# Get the distance to this feature (it is not necessarily the nearest one)
mindistance = pt.distance(nnfeature.geometry())

The above line is where I get an error. Python or QGIS can't find the distance between a point a line geometry. If nnfeature was a point geometry, it would work.

px = pt.x()
py = pt.y()
# Get all the features that may be closer to the point than nnfeature
closefeatureids = spacialIndex.intersects(QgsRectangle(px - mindistance,
                   py - mindistance, px + mindistance, py + mindistance))
for closefeatureid in closefeatureids:
    closefeature = featureID_Layer.getFeatures(QgsFeatureRequest(closefeatureid)).next()
    thisdistance = pt.distance(closefeature.geometry())

    # Python or QGIS doesn't like the line above either

    if thisdistance < mindistance:
        mindistance = thisdistance
        nnfeature = closefeature
        if mindistance == 0:
             break
#End of New NN logic

genericFeat2 = QgsFeature()
#convert from new NN logic
genericFeat2=nnfeature
idx = featureID_Layer.fieldNameIndex(fieldname)
closestFeatureFieldNameString=str(genericFeat2.attributes()[idx])
return closestFeatureFieldNameString


texts: Where can I find Population data for all UK Counties?

texts: Where can I find population data (male, female, total) for all UK counties for last 7 decades, i.e. 1950,1960,1970,1980,1990,2000,2010. I tried my luck with Office for National Statistics and UK National Statistics, but all in vain.


texts:
I'm starting to mess about with the solar radiation packages of GRASS and ArcGIS and have become rather stuck on the GRASS side of matters. I have a DSM area, and that is feeding into GRASS well enough, however when processing on the default values it falls over, not outputting anything. I'm using the command...

g.region rast=elevation.dem
r.sun -s elev=elevation.dem day=180
r.colors col=grey
d.his i_map=beam_map h_map=elevation.dem

In the R.Sun man page, the following command is recommended:

g.region rast=elevation.dem
r.sun -s elev=elevation.dem slop=slope asp=aspect beam=beam_map day=180
r.colors beam_map col=grey
d.his i_map=beam_map h_map=elevation.dem

However here, where am I supposed to get the slope, aspect and beam from?


texts: Split Parts Tool Error: Geometry invalid

texts:
I have trouble splitting a multipart feature and can't figure out what I'm doing wrong.
I have a multipart feature of Islands which all belong to the same country, so they should stay a multipart feature. However, in some cases the shapefile wasn't digitized properly so there is one larger polygon that in reality should be a few little ones. I tried using the "Split Parts" tool and it works for the first split so that I have two instead of one polygon, but when I try to split it again I get an error saying the geometry is invalid and I need to fix it before I can split.


texts: Batch append features with custom attributes

texts:
I have a directory with 2,000+ subdirectories that contain shapefiles submitted by project managers through a web form. Some contain additional subdirectories that hold the shapefiles as follows:

L:\uploads\102154682654842\14-52-63.shp or \102151486584621\project\14-59-22.shp

The submitted shapefiles have no attribute field standards, so the field names are comletely random.

I would like to append the shapefiles to a feature class and assign the directory folder name (which is a unique number) to the master feature class ID field and the shapefile name project.shp to the master shapefile name field, such that the master feature class attribute table looks like this:

        ID       |   Name   |
 ----------------|----------|
 102154682654842 | 14-52-63 |
 ----------------|----------|
 102151486584621 | 4-59-22  |
 ----------------|----------|

I can parse those strings with os.walk

And I can append the shapefiles to the master feature class using something like this:

import arcpy, os

inputFolder = r"L:\Uploads"
#filelist = list()
errorcount = 0
passcount = 0

# traverse root directory, and list directories as dirs and files as files
for root, dirs, files in os.walk(inputFolder):
    path = root.split('/')
    arcpy.env.workspace = root
    features = arcpy.arcpy.ListFeatureClasses(feature_type='polygon')
    #filelist.append(features)

    try:
        passcount += 1
        arcpy.Append_management(features, "L:\Test.gdb\Polygon", "NO_TEST","","")
#### ASSIGN path[index] to featureclass.Id and feature to featureclass.Name HERE! ####
        print "Append successful: ", path, features
        print "successful count: ",passcount

    except Exception, e:
    errorcount += 1
    print "Failed to append: ", path, features
    print 'for the following reasons: '+ str(e)
    print "error number: ",errorcount

Is there a way to add these attributes to the master feature class, through field mapping or other data management methods via arcpy?


texts: What is the single most influential book every GIS professional should read?

texts:
If you could go back in time and tell yourself to read a specific book/text/guide/tutorial at the beginning of your career as a gis-developer, which book would it be?

I expect this list to be varied and to cover a wide range of things.

This same question was once asked at stackoverflow for programming topics, and it has been considered by some as "one of the most important question ever asked on stackoverflow".
So I decided to bring this question to serve our needs here at gis.stackexchange.


texts: ArcGIS JS - Feature Layer selection symbol before layers loaded

texts:
I have a feature layer that is dynamic and could pull a

geometryType: esriGeometryPoint

or

geometryType: esriGeomertyPolygon

I was wondering how would I know when to set the selection symbol to

new esri.symbol.SimpleMarkerSymbol(...)

or

new esri.symbol.SimpleFillSymbol(...)

Before I was able to just set the selection symbol to new esri.symbol.SimpleMarkerSymbol(...) and it would fill the polygon symbol.  But now something seems to have changed where it now just places a circle or marker symbol in its place.

I've tried to set the selection symbol on featureLayer.on("load", function (evt) { ... }), but it just seems to use the default color of the layer.

Thanks in advance.


texts:
Is it possible in osm2po to additionally consider surfaces as edge weights along with the highways?

I've tried the following:

wtr.tag.highway.pedestrian = 1,1,5, foot

wtr.tag.highway.footway = 1,1,5, foot

wtr.tag.surface.paving_stones = 2,1,1, foot

wtr.tag.surface.asphalt = 2,1,4, foot


texts: Is there a way to force leader lines to polyline centroids in ArcGIS Desktop 10?

texts:
I'm making a pretty standardized set of maps that include alphabetical labels for specific stream segments.  As it stands, all my label leader lines snap to the downstream end of my polylines.  My current workaround is to convert labels to annotation and hand-place them but this is getting to be a real time sink.  Thus I ask: is there a setting I'm missing, or other way to get leader lines to snap to a line centroid?

Full disclosure: [edited to add] Running ArcEditor license level.  I don't have the Maplex extension and am not likely to get it.


texts: How to get dynamic county information into mapbook?

texts:
I have access to ArcMap 10.0 (ArcInfo/Workstation) on Windows 7 SP1.

My usual scenario is to produce a map book having 3 shapes:

1 polyline (a road track about 200-300 km without any useful attribute data)
1 polygon (boundaries with all the county data like county names, numbers etc. in the attribute table)
1 polygon (road track based Strip Map Index Feature (SMIF) to get the DDP running with attributes generated by ArcMap)

All works fine and ArcMap produces around 130-150 ISO A4-PDF pages per map book at a given scale.

But how to get the county names fieldbased on each single page without some manual action?

I know that ArcMap only uses the SMIF-Layer with fixed attribute data as a valid source for DDP, so no Joins/Relates to the attribute table are allowed. And adding new fields manually to SMIF table and filling them actually works, I did that last time, but, you know, it stinks somehow.

So, how to get the attributes from boundaries.shp into SMIF.shp to provide dynamic data for each page in DDP?


texts: Relative positional accuracy vs. registration error - what is the correct terminology?

texts:
I have several aerial orthorectified images that I am using for vector-based (points) change detection. I have calculated RMSE using sets of points on static landscape features in the images (e.g., an RMSE value for images from flight 1 and flight 2). In this case, absolute accuracy is not important, or at least much less important than relative accuracy.

What is the appropriate terminology to use when reporting these values?
Relative (horizontal) positional accuracy?
Co-registration error?
What is the difference here?

I've always assumed that registration error implies that one image was registered to another, but after reading some literature where this is unclear now I am not sure.


texts: How to retrieve the mapPoint behind a PictureMarkerSymbol?

texts:
I am using ArcGis Javascript API 3.2. After instantiating a Point, PictureMarkerSymbol, attribite and an InfoTemplate, I use this code to plot a point on the map.

var graphic= new esri.Graphic(candidatePoint,candidateSymbol,attr,infoTemplate);
map.graphics.add(graphic);

I want to retrieve the point behind the graphic when it is being clicked. How will I do it?

Thanks,
Prasanna


texts: Add GeoJSON features from an API call in OpenLayers 4

texts:
I'm trying to create a little demo to show the data from our API into a map but I have all the time the error:

Uncaught TypeError: feature.getId is not a function

Here my code, is just basic for the moment, we just need to display the data:

var map = new ol.Map({
  target: 'map',
  layers: [
    new ol.layer.Tile({
      source: new ol.source.OSM()
    })
  ],
  view: new ol.View({
    center: ol.proj.transform([-0.92, 52.96], 'EPSG:4326', 'EPSG:3857'),
    zoom: 8
  })
})

var appSource = new ol.source.Vector({
  projection: 'EPSG:27700',
  format: new ol.format.GeoJSON()
})

var appLayer = new ol.layer.Vector({
  source: appSource,
  style: new ol.style.Style({
    fill: new ol.style.Fill({
      color: 'rgba(255, 255, 255, 0.1)'
    }),
    stroke: new ol.style.Stroke({
      color: 'rgb(85, 172, 238)',
      width: 2
    })
  })
})

$.ajax({
  type: 'GET',
  url: 'https://www.traffwebdev.uk/tro/v1/restrictions?epsg=&x=383666+&y=398197&area=&radius=20&output=',
  dataType: 'json',
  success: function (data) {
    // reads and converts GeoJSon to Feature Object
    var features = (new ol.format.GeoJSON()).readFeatures(data.geometry)
    appSource.addFeature(features)
  }
}).then(() => map.addLayer(appLayer))

You can use this URL to see the returned data from the API, supposedly the ID is already there and in the correct place

https://www.traffwebdev.uk/tro/v1/restrictions?epsg=&x=383666+&y=398197&area=&radius=20&output=


texts:
I'm making the transition from GIS to PostGIS and am stuck on this spatial query.

I have background geometry (purple, outlined red) which has an overlay of islands (blue, below); I need the area of only the islands within the red line. In GIS this would be a clip, using the purple shape as my mask.

So I query using ST_Within() and obtain the green islands, which is great, but - logically - I'm missing those which are not fully within the mask:

And I can still find those islands overlapping the boundary of my polygon using ST_Overlap() to find these pink polygons:

But is there an efficient way to get the area of the islands only within the red outline? This will cut a percentage of the island polygons, I understand, those will be associated with the adjacent background polygon.


texts:
Boston housing data provided by GeoDa center is a well-known dataset.
GeoDa Center provides this data in shapefile format. A metadata file included in the download indicates that the data use projection UTM 19N. However the shapefile does not have a projection defined (no accompanying .prj file), and I am unable to load it correctly.

I would like to overlap the Boston data with the census tract data available at Census Bureau using ArcGIS for Desktop.

Do you know how I can do this?


texts: How to get the data on screen after movened event is triggered

texts:
I am displaying 100 cities on the map.

After zoom or move the map, only e.g. 9 items are displaying on the current screen.

So how can I get these data from my dataset? Especially the object id?

Following is the code to create layer:

var sql = new cartodb.SQL({ user: 'achilles' });
        var query = "SELECT * FROM mi_munv5 where objectid in (" + ids + ")";
        layer.setSQL(query);


texts: Creating Story Map based on image keywords?

texts:
I would like to create an interactive map that displays photos based in their geotags but could also be sorted by other metadata components like keywords. The keywords are mainly based on species, so ideally the user can click on the species and view their distribution within the map.

I think I am looking for a layout similar to Esri's Crowdsource Story Map .

If not Story Map, I am open to other mapping platforms.


texts: How to fix encoding problems with Turkish characters?

texts:
Possible Duplicate:
  How to read Greek fonts (ISO-8859-7) in shapefile attributes within QGIS 1.8.0?

My shapefile contains turkish (latin) letters. (see: http://www.5goldig.de/Tuerkische_Tastatur/Tuerkisch-Buchstaben.html)
The encoding of the shapefile is UTF-8 and I also tried ISO-8859 and several others but the letters are still like that:

ApÄ±Åkar GeÃ§it -> Apişkar Geçit

Is it because of the encoding or is there a other reason?


texts: How do you create a new GRASS mapset in QGIS 2.12? The GRASS panel has only the Open GRASS Tools button.I don't see any way to actually create a GRASS mapset in QGIS.


texts:
This is my script where I am trying to use selectby_attribute
within a model, based on a field, F*.  There are 6 fields, so I want 6 derived layers.  Then I am going to use bearing and distance tool on each of these output layers.  I cannot pass query to the model.  How do I call a script from my 1st model, then call a model from this script?

import arcpy, sys, string, os

#arcpy.env.workspace = arcpy.GetParameterAsText(0)
in_FC = arcpy.GetParameterAsText(0)
#FList = arcpy.GetparameterAsText(1)
#F_list = FList.split(";")
#arcpy.ListFields(in_FC, "Rev_2016_F*")

#arcpy.AddMessage(tableName)
arcpy.AddMessage("The table %s exists" % (in_FC))
CurrentWorkspace = r"D:\Fort_Belvoir\Conservation\Conservation.gdb"

#in_Table = "D:\Fort_Belvoir\Conservation\Conservation.gdb\QryTab4_CopyRows"

#tbx = arcpy.AddToolbox("C:/Temp/MyToolbox.tbx")
tbx = arcpy.AddToolbox(r"D:\Fort_Belvoir\Conservation\PIF.tbx")
#arcpy.AddMessage("this is toolbox %s" % (tbx))

# Get a list of field objects
#
num = 0
fields = arcpy.ListFields(in_FC, "Rev_2016_F*", "LONG")

for field in fields:

    # Check the field name, perform a calculation when finding the field 'Flag'
    #
    # Local variables
    anum = str(num)
    thetab = "F" + anum
    arcpy.AddMessage(thetab)
    #arcpy.SetParameter(1, SR)
    theFld = field.name
    arcpy.AddMessage(field.name)
    arcpy.AddMessage(theFld)
    arcpy.AddMessage(str(num))
    anum = str(num)
    arcpy.AddMessage(anum)
    thetab = "F" + anum
    arcpy.AddMessage(thetab)

    out_FC = r"D:\Fort_Belvoir\Conservation\Conservation.gdb\PIF" + os.sep + thetab

    try:
        arcpy.AddMessage(out_FC)
        arcpy.SetParameter(1, theFld)
        arcpy.SetParameter(2, out_FC)
        tbx.Model2_Toolbox(theFld, out_FC)

    except:
    arcpy.AddMessage("Failed")

    arcpy.AddMessage(out_FC)
    num = num + 1


texts: How to Turn Off SFCGAL Geometry Validation?

texts:
SFCGAL documentation (https://oslandia.github.io/SFCGAL/documentation.html)
states that all geometry processing algorithms by default check the validity of input geometries, but also come in a second variant "where no actual check is made (taking a NoValidityCheck tag as additional parameter)".

However, I can't find any examples or detailed documentation for queries that use a NoValidityCheck tag.

What is the syntax for this and does it actually work with the SFCGAL extension to postgresql?


texts: How to customize grids in Map Composition in QGIS?

texts: Please somebody help me out with Map Composition. In map composer how to put in geographical grid (not UTM grid)? And how to put only frame with ticks and annotation (not entire grid across a map)?


texts: ST_Intersection between MultiPolygon tables leads to gaps, missing intersections

texts:
I am trying to intersect two MultiPolygon tables, which are in the same UTM projection, and calculate the area of their intersections. I am getting some surprising results when I examine the intersections, which I retain in a new table. Specifically, each table contains either 2000 or 2010 Census block group boundaries. Here, I find the intersections with a WHERE clause, calculate the weight associated with the area of that intersection, and put the results into a new table.

SELECT c2000.bkgpidfp00 AS fips2000, c2010.geoid10 AS fips2010,
  ST_Area(ST_Intersection(c2000.geom, c2010.geom)) / ST_Area(c2000.geom) AS weight,
  ST_Intersection(c2010.geom, c2000.geom) AS geom
  INTO xwalk
  FROM census2000_block_groups c2000, census2010_block_groups c2010
 WHERE ST_Intersects(c2000.geom, c2010.geom) AND NOT ST_Touches(c2000.geom, c2010.geom);

Problem is, there are gaps in the final result. Not all of the intersections are present! See images below.

Here is the intersection--gaps are in white.

Here are the 2000 (blue) and 2010 (red) block group boundaries--you can see there are other intersections that should have been found!

Finally, just to prove that these intersections make sense, I used QGIS to calculate the intersections as well--and it worked! **Below is the correct intersection from QGIS ("Vector" > "Geoprocessing" > "Intersect").

Any ideas why my PostGIS ST_Intersection example is failing to produce all of the possible intersections?

In the interest of reproducibility, here's a shell script to download and subset the data from the U.S. Census:

wget ftp://ftp2.census.gov/geo/pvs/tiger2010st/26_Michigan/26/tl_2010_26_bg00.zip
wget ftp://ftp2.census.gov/geo/pvs/tiger2010st/26_Michigan/26/tl_2010_26_bg10.zip

unzip ~/Downloads/tl_2010_26_bg00.zip
ogr2ogr -f "ESRI Shapefile" -where "COUNTYFP00 IN ('099', '125', '163')" -s_srs "EPSG:4269" -t_srs "EPSG:32617" census_block_groups_2000_utm17n.shp tl_2010_26_bg00.shp

unzip ~/Downloads/tl_2010_26_bg10.zip -d $WD
ogr2ogr -f "ESRI Shapefile" -where "COUNTYFP10 IN ('099', '125', '163')" -s_srs "EPSG:4269" -t_srs "EPSG:32617" $WD/census_block_groups_2010_utm17n.shp tl_2010_26_bg10.shp

To insert into PostGIS:

shp2pgsql census_block_groups_2000_utm17n.shp public.census2000_block_groups | psql -h localhost -d detroit
shp2pgsql census_block_groups_2010_utm17n.shp public.census2010_block_groups | psql -h localhost -d detroit

Software versions are:

QGIS: 2.8.1-Wien
GDAL/OGR: 1.11.2
GEOS: 3.4.2-CAPI-1.8.2
Proj.4: 480
PostgreSQL: 9.3.6
PostGIS: 2.1.0


texts: ogr2ogr spatial query doesn't work in process

texts:
I tried exporting data selecting it from postgres with spatial query. The argument does work within osGeo4W Shell but I'm trying to use it from my C# console app. Everything seems too work until I add the spatial part. I have no idea what could be the issue here, but I'm still a beginner. Can someone take a look and help me figure out what went wrong here?
// assemble input
string output = "C:\\Users\\martijn\\Documents\\output.gpkg";
string connectionPostgis = " PG:\"dbname = db_naam host = localhost user = postgres password = password port = 5432 \" ";
string sqlRequest = "-sql \"SELECT naam, geom,ST_Distance(ST_GeographyFromText(\'POINT(4.9235515 52.3770089)\'),geom) FROM public.table \"";
string argument ="ogr2ogr -f \"GPKG\" " + output + connectionPostgis + sqlRequest;

// process
Process process = new Process();
process.StartInfo.FileName = "C:\\OSGeo4W64\\OSGeo4W.bat";
process.StartInfo.Arguments = argument;

process.Start();


texts: unable to install qgis 1.9 on Lubuntu 13.04

texts:
I tried to install qgis 1.9 in my laptop. I use Lubuntu 13.04 (64bit).
I used the repository

deb     http://qgis.org/debian-nightly raring main
deb-src http://qgis.org/debian-nightly raring main

and I didn't find any error during the installation.

When I try to start qgis I obtain this message:

  /usr/bin/qgis.bin: error while loading shared libraries: libgdal1.7.0.so.1: cannot open shared object file: No such file or directory

If I try

sudo apt-get update
Sudo apt-get upgrade

I obtain those messages:

Configurazione di qgis-providers-common (1.9.0+git20130529+e0406e7~raring1)...

/usr/lib/qgis/crssync: error while loading shared libraries: libgdal1.7.0.so.1: cannot open shared object file: No such file or directory
dpkg: errore nell'elaborare qgis-providers-common (--configure):
 il sottoprocesso installato script di post-installation ha restituito lo stato di errore 127
dpkg: problemi con le dipendenze impediscono la configurazione di qgis-providers:
 qgis-providers dipende da qgis-providers-common (= 1.9.0+git20130529+e0406e7~raring1); comunque:
  Il pacchetto qgis-providers-common non è ancora configurato.

dpkg: errore nell'elaborare qgis-providers (--configure):
 problemi con le dipendenze - lasciato non configurato
dpkg: problemi con le dipendenze impediscono la configurazione di qgis:
 qgis dipende da qgis-providers (= 1.9.0+git20130529+e0406e7~raring1); comunque:
  Il pacchetto qgis-providers non è ancora configurato.

dpkg: errore nell'elaborare qgis (--configure):
 problemi con le dipendenze - lasciato non configurato
dpkg: problemi con le dipendenze impediscono la configurazione di qgis-plugin-globe:
qgis-plugin-globe dipende da qgis (= 1.9.0+git20130529+e0406e7~raring1); comunque:
Segnalazione apport non scritta poiché il messaggio di errore indica la presenza di un fallimento precedente.

Segnalazione apport non scritta poiché è stato raggiunto il valore massimo di MaxReports

  Il pacchetto qgis non è ancora configurato.

dpkg: errore nell'elaborare qgis-plugin-globe (--configure):
 problemi con le dipendenze - lasciato non configurato
dpkg: problemi con le dipendenze impediscono la configurazione di qgis-plugin-grass:
 qgis-plugin-grass dipende da qgis (= 1.9.0+git20130529+e0406e7~raring1); comunque:
  Il pacchetto qgis non è ancora configurato.

dpkg: errore nell'elaborare qgis-plugin-grass (--configure):
 problemi con le dipendenze - lasciato non configurato
Si sono verificati degli errori nell'elaborazione:
 qgis-providers-common
 qgis-providers
 qgis
 qgis-plugin-globe
 qgis-plugin-grass
E: Sub-process /usr/bin/dpkg returned an error code (1)

I apologize for the Italian messages.

Could you help me?

Thank you


texts: extract separate nodes in Postgis equivalent QGIS

texts:
I want to extract nodes from polygon in Postgis like it QGIS "extract nodes" tool does. I've red the answered question about this in What's the PostGIS equivalent to Extract Nodes in fTools and used ST_DumpPoints(geom) function. But the result contains both the geometry and a path telling where in the polygon the point belongs. (ST_DumpPoints(the_geom)).geom gives the geometry (points), OK.
The problem is that I can't address to single node from geometries, because tool ST_DumpPoints return integer indexes (path) of the polygon the point belongs to. In short, when I visualize my dump points in QGIS, I have no possibility to select one single node - all nodes of polygon have been filtered (picture 2), but selected is only one. As a consequence, I can't implement any spatial analysis with every separate node (e.g. verification about other layer contains any of this nodes).

How I can simply get the points without connection to the polygons?


texts: Find the right values for bounds and origin for Proj4Leaflet

texts:
I'm using Proj4Leaflet to use tile layers with a different CRS (EPSG:25832) in my leaflet project.

Now I am somehow lost in finding the correct values for the bounds and origin of the tiles I want to use (https://sgx.geodatenzentrum.de/wmts_topplus_open/tile/1.0.0/web/default/EU_EPSG_25832_TOPPLUS/{z}/{y}/{x}.png). I've constructed this tile request by inspecting the capabilities of this WMTS (here is a link to the metadata of the service).

Thanks to this post, I've already figured out how to calculate the values for the resolution.

var epsg25832 = new L.Proj.CRS(
  'EPSG:25832',
  '+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs', {
    origin: ? ? ? ,
    bounds: ? ? ? ,
    resolutions: [
      4891.969810252,
      2445.984905126,
      1222.9924525616,
      611.4962262808,
      305.7481131404,
      152.87405657048,
      76.43702828524,
      38.21851414248,
      19.109257071296,
      9.554628535648,
      4.777314267824,
      2.388657133912,
      1.194328566956,
      0.597164283478
    ],
  }
);

But I have no clue how to find the right values for bounds and origin. I can only find <ows:WGS84BoundingBox> in the capabilities of the WMTS. Does anyone know how to get the right values?


texts: How to change the font size of the color bar of a GeoPandas choropleth plot

texts:
When I try to use the legend_kwds argument to change the font size of my colorbar, I keep getting this error:

TypeError: init() got an unexpected keyword argument 'fontsize'

ax = df.plot(figsize=(20,16), alpha=0.8, column='value', legend=True, cmap='OrRd', legend_kwds={'fontsize':20})

plt.show()

Does anyone know how I can increase the font size of the colorbar with GeoPandas? I can't seem to find a keyword that works. I'm using GeoPandas 0.8.1 and Matplotlib 3.3.1.


texts: How many stereo satellite imagery required to prepare DEM of an area?

texts: I want to know how to find the nos. of scene of stereo CARTISAT-1 data required to prepare DEM of an area(say 84000 sq km)? any suggestion/solution/calculation?


texts: Point “inside” many polygons (leafletjs, turfjs, geojson)

texts:
I'm trying to use "inside" function of turf js and leafletjs. I've prepered the basic example with one geojson point and one polygon but I cannot do the same with more than two polygons (geojson). Any ideas?

<!DOCTYPE html>

    <title>Mobile tutorial - Leaflet</title>

    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0">

    <link rel="shortcut icon" type="image/x-icon" href="docs/images/favicon.ico" />

    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.0.1/dist/leaflet.css" />


texts: How to automate point-in-polygon queries in ArcView?

texts:
Layer1 = point shapefile (multiple points, geocoded)

Layer2 = polygon (1 attribute; (1,2,3...))

I want a non-gis person to be able input a single address and be told the attribute of the polygon which it falls within.

What is the simplest way to accomplish this?

The user will be doing this in arcview.  I use ArcInfo.


texts:
Currently running Version 2.14.13

I'm having trouble getting base maps to appear when I go to print a map in QGIS Cloud. I've been using Google Hybrid from the OpenLayers plugin, and while it appears when I pull up my web map, it disappears when I attempt to use the print feature. I realized OpenLayers is outdated so I tried to use the QuickMapServices plugin instead.

I've tried using OSM Standard from the QuickMapServices plugin. However, when I go to publish my map in QGIS, I get an error message stating "Unsupported Layers Found" in reference to the OSM Standard layer. I know there is a way to add a service within the QuickMap setting, but this is a process that I'm not familiar with, nor do I know if this is something that will even make a difference here.
I realize I could just upload some aerial photos, but the size of those files doesn't really make it feasible with the amount of DB storage I have left. Any ideas of a solution?


texts: I have a single layer of polygons (drawn around every zip code). I have a CSV file with territory numbers and the zip codes in each territory. I joined these and was able to color each territory (multiple zip codes) instead of each individual zipcode. What I'd like to do now is get rid of the zip code boundaries to only have the territory boundaries, but I can't figure out how to do this.


texts: Converting KMZ file with png images to shapefile?

texts:
I have a KMZ file with data cell phone antenna coverage.

I try to convert with tool kml to layer but it converts only points. That information with coverage disappeared. I tried to use ET but it also returns only points.

Any ideas?


texts: GeoPandas 0.6.1 installed instead of 0.7.0 in conda windows 10

texts:
I have problems using the .clip utility from the GeoPandas package. I found in this post that the problem might be due to the GeoPandas version, since this utility is only available in > 0.7 version.

https://stackoverflow.com/questions/61505660/geopandas-has-no-attribute-clip

I tried conda install geopandas, but always I receave the notification that the version 0.6.1 will be installed, and not the version 0.7. Is there a way to do it through conda and not pip?


texts:
I am trying to draw ellipse on the client side using openlayers with some mathematical functions from drawing ellipse with svg. I have achieved drawing it but not rotation what I have written for parameter and I am not sure about axis length too.

When I measure semi-major or semi-minor axis length, I am not getting the true length what I have given for 4.8 * 1000 meter.

I know PostGIS 2.1 has ST_MakeEllipse function for this but I want it on the client side.

I am sharing javascript code, reference image and result image with you. Any help would be appreciated.

Reference Image:

Codes:

var radiusX= 4,8 * 2 // Semi-major axis length - km
var radiusY= 2,5 * 2 // Semi-minor axis length -km
var rot = 125 //Rotation

var xCent = 14.7495;
var yCent = 47.4640;

var z = [];
var list = [];

for (var i = 0 * Math.PI; i < 2 * Math.PI; i += 0.01 ) {
    xPos = xCent - (radiusX* Math.sin(i)) * Math.sin(rot * Math.PI) + (radiusY* Math.cos(i)) * Math.cos(rot * Math.PI);
    yPos = yCent + (radiusY* Math.cos(i)) * Math.sin(rot * Math.PI) + (radiusX* Math.sin(i)) * Math.cos(rot * Math.PI);

    if(i == 0 ) {
        z.push( new OpenLayers.Geometry.Point(xPos, yPos).transform(map.displayProjection, map.projection));
    } else {
        list.push(new OpenLayers.Geometry.Point(xPos, yPos).transform(map.displayProjection, map.projection));
    }
};

var linear_ring = new OpenLayers.Geometry.LinearRing(list);
var polygonFeature = new OpenLayers.Feature.Vector(new OpenLayers.Geometry.Polygon([linear_ring]), null, null);
vector.addFeatures([polygonFeature ])

Result:


texts: pour point user click input over map for watershed delineation

texts:
I have a drainage density map and i want to interactively retrieve watershed by clicking a point on map. Once I zoomed to an area of interest and click the highest accumulation point of my choice then i get in return watershed only encompassing that outlet. I am wondering if anyone had already done that manually or automatically. Any suggestion to do it will be helpful.

Thanks.


texts:
I have begun learning how to manipulate LAS data in python and wanted to see how others handle LAS files.  I would like to read the points (I am using a numpy array), and filter out classes 1 and 2 (unclassified and ground) to a separate array.  I have the following code but cannot seem to get the points filtered.

# Import modules
from liblas import file
import numpy as np

if __name__=="__main__":
    '''Read LAS file and create an array to hold X, Y, Z values'''
    # Get file
    las_file = r"E:\Testing\ground_filtered.las"
    # Read file
    f = file.File(las_file, mode='r')
    # Get number of points from header
    num_points = int(f.__len__())
    # Create empty numpy array
    PointsXYZIC = np.empty(shape=(num_points, 5))
    # Load all LAS points into numpy array
    counter = 0
    for p in f:
        newrow = [p.x, p.y, p.z, p.intensity, p.classification]
        PointsXYZIC[counter] = newrow
        counter += 1

I have seen arcpy.da.featureClassToNumpyArray, but I did not want to import arcpy nor have to convert to shapefile.

How else can I filter/read LAS data into a numpy array?


texts: Shapefile partly hidden in QGIS but fully visible in Geosecma

texts:
I'm working in QGIS Hannover and have added some vector shapefiles to my project. Some of the data in the files won't show up in the viewport, leaving "holes" in the map. The same file contains both visible and "invisible" data.
A colleague of mine opened them in Geosecma (an ArcMap module) for me. She discovered that the data not showing up was displayed from empty categories in the attribute table. When she changed this, she could view the complete dataset in the viewport. So there's nothing wrong with the data.
I went back into QGIS and tried displaying from every available attribute, using Symbology. I even tried using the exact same file as my colleague instead of my own copy. It did not solve the problem.
I'm at a loss. Any ideas?


texts:
After read this How merge osm files with same ID objects? I can't get any conclusion so please don't flag as duplicate.

I have some PBF files came from NASA SRTM ( http://katze.tfiu.de/projects/phyghtmap/phyghtmap.1.html )

My problem is phyghtmap have its own ID numbers and I believe start from same number in every PBF file.

I need to import these PBF files to my database (appart from OSM tables so official OSM ID is not a problem) using osm2pgsql.

For each HGT file I have a PBF and its own ID range that will conflict with the next files and override previous data in database.

How can I import all these files incrementing the ID for every file?

Tools I have: osmosis, osm2pgsql, osmium and phyghtmap.


texts: Why Spatial autocorrelation (Global Moran's I) is producing a p value of 0?

texts: I am undertaking spatial autocorrelation (Global Moran's I) in ArcMap 10.3. The p value returned is 0, does this just mean its extremely significant or is this an incorrect result?


texts: Date and geometry from JSON objects to Google Earth Engine?

texts:
I would like to know if there is any way to get the geometry and date data from a JSON object. I'm putting the coordinates and dates by hand but all my data comes from JSON files that have either the date and coordinates (https://0162ago2m5.execute-api.us-east-2.amazonaws.com/prod/life/133/)

Currently i'm using this code :

var geometry = ee.Geometry.Point([67.42467260724823, 21.057842149424154]);

var rrs = ee.ImageCollection("NASA/OCEANDATA/MODIS-Aqua/L3SMI")
  .filterDate('2018-1-04', '2018-1-05')
  .select('Rrs_488', 'Rrs_531', 'Rrs_547', 'Rrs_555')
  .filterBounds(geometry);

rrs = rrs.toList(rrs.size());

var reduceR = function(image) {
  image = ee.Image(image);
  var reducedVal = image.reduceRegion({
    reducer: ee.Reducer.first(),
    geometry: geometry,
  })
  return ee.Feature(geometry).set(reducedVal);
};

var mRrs=rrs.map(reduceR);
mRrs = ee.FeatureCollection(mRrs);

Export.table.toDrive({
  collection: mRrs,
  description: 'rrs',
  fileNamePrefix: 'rrs',
  fileFormat: 'CSV',
  selectors: ['Rrs_488', 'Rrs_531', 'Rrs_547', 'Rrs_555']
});


texts: How to make points with fixed distance and orientation from another points layer?

texts:
I have a points layer with hundreds of features at variable distance from each other. Now what I want to do is create 4 new points for each feature of this layer at a fixed distance (30m) AND orientation (exactly North, East, South, West).

Something like this (forgive my draw skills):

Where grey points are from input and the color ones are the output.

I was willing to do in QGIS, but can be also PostGIS or another solution.


texts: Getting a frequency matrix of overlaying raster values

texts:
I have two raster datasets (a World Soil Map and Global Horizontal Irradiance), which I merged together in QGIS in order to get values of both channels for every pixel.

The Soil Map has 30 soil type classes and the GHI raster has been restricted to 16 classes.

I would like to get an overview of the frequency of each combination (30x16).

The next example shows the values of a random point. As the upper table is simplified, here you get an actual point of the raster: Channel 1 shows the soil type (57 represents Gleysols) and channel 2 gives the value of irradiance (Long-term average of daily sum <1752 kWh/m^2):

Do you have any recommendations how to process a frequency matrix of this?
I also want to do this using vector data and compare the results.
(My R and Python-Skills are rather limited, but I'm willing to learn)


texts: What is the correct way to query for traffic_signs in Overpass in an area?

texts:
I am new to all this stuff.
I currently want to get all the traffic signs in an area. I am running this query in overpass turbo, but nothing is returning back to me.
/*
This is an example Overpass query.
Try it out by pressing the Run button above!
You can find more examples with the Load tool.
*/
node
  [traffic_signs='*']
  ({{bbox}});
out;


texts:
Having trouble here: I'm trying to use python to write selective label expressions. I have several concentrations and depending on its "_Ex" number I either want it labeled red, not red, or not labeled at all. I have started with just one trying to set it to red. when i put this into the lable expression it says line one Findlabel not defined.

Any help would be appreciated. -Acrain  (Indentations are not showing up correctly here, but I don't think thats my problem)

def FindLabel ([cis12DCE] , [cis12DCE_Ex]):

if short ( [cis12DCE_Ex] ) = 1:
    return "<clr red = "255">"+"cis-1,2-Dichloroethylene  "+[cis12DCE]+"</clr>
else:
    return "cis-1,2-Dichloroethylene  "+[cis12DCE]


texts: Calculating midpoint for snapping using OpenLayers?

texts:
How to calculate midpoint of all lines in layers, and add to it layer for midpoint snapping, using OpenLayers 2.13?

I've tried the following code, but does not work:

var features=vectors.features;
var len=vectors.features.length;
for (var i=0 ; i &lt; len ; i++){
    for(var j=0 ; j&lt;features[i].geometry.getVertices().length ; j++){
        points.push(features[i].geometry.getVertices()[j]);
    }
}
i=0;
while(i&lt;points.length){
    xA=points[i].x;
    yA=points[i].y;
    xB=points[i+1].x;
    yB=points[i+1].y;
    xC=(xA+xB)/2;
    yC=(yA+yB)/2;
    midPoints.push(new OpenLayers.Feature.Vector(new OpenLayers.Geometry.Point(xC,yC)));
    i++;
}


texts: Symbolizing raster layer with Python QGIS

texts:
I have a simple Python code that makes some calculations on rasters and give me some output rasters; I have also a style saved to color rasters based on pixel values.
Is there an option to load the output (single-band) rasters with the color style I want in stead to have them in black and white?
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
b=[12.158, 13.094, 15.961, 19.553, 24.669, 29.16, 31.805, 31.615, 27.284, 22.286, 17.026, 13.31]
s=['Tmaxjan', 'Tmaxfeb','Tmaxmar','Tmaxapr','Tmaxmay','Tmaxjun','Tmaxjul','Tmaxaug','Tmaxsep','Tmaxoct','Tmaxnov','Tmaxdec']

path = '/Users/macbook/Desktop/TESI/PROGETTO/raster_T_max/'

list=[0,1,2,3,4,5,6,7,8,9,10,11]
for i in list:

  rasterfile = QgsRasterLayer('/Users/macbook/Desktop/TESI/PROGETTO/DEM CAPITANATA.tif')

  entries = []
  ras1 = QgsRasterCalculatorEntry()
  ras1.ref = 'ras@1'
  ras1.raster = rasterfile
  ras1.bandNumber = 1
  entries.append( ras1 )

  new_path = path + s[i]+ '.tif'

  temperature = str(b[i]) + "+" + '(-0.007)' + '* ras@1 '

  Temp = QgsRasterCalculator(temperature, new_path, 'GTiff', rasterfile.extent(), rasterfile.width(), rasterfile.height(), entries )

  Temp.processCalculation()
  iface.addRasterLayer(new_path)


texts: Set property disappears in ImageCollection mapped function (Earth Engine)

texts:
I am creating a ML classifier to detect floods, but I need to find imagery to train it with. To do this, I am using a Fusion Table with thousands of flood events, their began/ended dates, and polygons of their extent.

I want to quickly find what imagery is available for these floods, then see if it is suitable for training. In GEE, I created a function to map over the fusion table, find Landsat image collections for the date range and bounds of a flood event, then put them together. The "putting together" part is causing me trouble.

I want to be able to set a flood ID property ('dfoID') for the images/image collections. That way I can see which flood they belong to later. Ideally, all the images from each flood event would be grouped into one image collection. However, because I map a Feature Collection, the output is a Feature Collection. I flatten this and then get one image collection but none of the dfoIDs remain on the images.

Here's the link to my code:
https://code.earthengine.google.com/3c31496721f4a01c4d9c77933e408266

Any ideas for 1. Yielding an image collection of image collections; Or 2. Making the dfoIDs stick to either the individual images, or to an image collection of images for each flood event?


texts: How to use variables in arcpy tools which quote the variable name

texts:
I'm having a frustrating time trying to use the Create TIN tool within Python, when the input featureclass name is a variable.

The help file entry for Create TIN shows a hard-coded path to the featureclass "points.shp" in the stand-alone example:

arcpy.CreateTin_3d("NewTIN", "....prj", "points.shp Shape.Z masspoints")

The script example uses {0}:

arcpy.ddd.CreateTin(outTin, sr, "{0} Shape.Z masspoints" .format(lasMP), "Delaunay")

It's unclear to me what the {0} represents (the first argument? - what if I don't have an argument, but only have access to a variable?).

The syntax section of the help file shows that the in_features can be specified as an array (of arrays?):

[[in_feature_class, height_field, SF_type, tag_value],...]

but I can't get this to work - it returns a generic "ERROR 999999: Error executing function".

So my question is, how can I use a variable as the featureclass name in the CreateTIN tool within Python?

Here is my non-functioning code:

fc = "path_to_featureClass"
tin = "path_to_output_TIN"
sr = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"

arcpy.CreateTin_3d(tin, sr, [fc + "fieldName Hard_Line <None>"], "DELAUNAY")

(where I have valid paths to fc and tin, and fieldName is a valid field on the layer.) I can successfully create a TIN using this dataset when using the tool manually.


texts:
I have a table in QGIS (using QGIS 2.4). I have created a form for adding attributes. When I add new objects on the layer, I also create it's attributes. Can I make attributes that I have added non-editable for users?

Maybe I can do it if I add a PostgresQL layer and create some restrictions?


texts: Import data, analyze spatially, and export list of points using Google Earth

texts:
I am trying to:

import data into Google Earth from an outside data source.
Then overlay a map and isolate which datapoints fall into that area.
I would like to then export a list of the data points to give to a person that works with our company.

Can I do this in Google Earth Pro?  Can you direct me to a resource that would walk me through this?  Is there a better solution for doing this type of work if it can't?


texts: Utilizing versioning with ArcGIS Collector

texts:
I have done a lot of research on the internet with this and I am unable to find a cut and dry answer. I am utilizing ArcMap 10.3.1. and PostgreSQL 9.3. I have been able to get everything set up between ArcMap and PostgreSQL. I created a enterprise geodatabase, imported some point data, registered the database with the server, and successfully published the data to ArcGIS Portal. Within Portal, I have made a map and developed an appropriate web application. With this I can successfully edit the data in the PostgreSQL database utilizing ArcGIS Collector.

I would like to utilize versioning in order to keep track of the changes that are in being made in the field. I have successfully created a version that is public and protected the DEFAULT version. When i do this, I am unable to sync the edits back to the PostgreSQL database.

Has anyone any insight or had any luck getting versioning to work with ArcGIS Collector?


texts: Internal problem using Google Earth Engine function

texts:
I need help solving this GEE internal error. I am trying to reduce my script lines in GEE. So, I tried to use a single function to compute polarization indexes.

However, the console shows a warning saying:

This is the script:

// Function to compute and add polarization indexes as bands
var addindexes = function(ima) {
  var NDPI = ima.select('VV').subtract(ima.select('VH'))
    .divide(ima.select('VV').add(ima.select('VH')));
  var NVHI = ima.select('VH').divide(ima.select('VV').add(ima.select('VH')));
  var NVVI = ima.select('VV').divide(ima.select('VV').add(ima.select('VH')));
  var VHrVV = ima.select('VH').divide(ima.select('VV'));
  var indexes = ee.List([NDPI,NVHI,NVVI,VHrVV]);
  return ima.addBands(indexes);
};

// Add index bands to image collection
var S1_indexes = allSeasons.map(addindexes);
print(S1_indexes);


texts:
Are Sentinel-2 data geometrically corrected or do they need to be georeferenced as an image to image?

What is the best software for processing?


texts: Nominatim - Run without reprocessing data

texts:
I have a Docker Nominatim instance up and running (https://registry.hub.docker.com/r/peterevans/nominatim/).
The command to start nominatim is: sudo -u nominatim /srv/nominatim/build/utils/setup.php --osm-file data.osm.pbf --all --threads 2
But that command reprocesses the whole OSM.PBF file every time I restart the docker container.
Is there a way to run nominatim without reprocessing everything again?
Something like "sudo -u nominatim /srv/nominatim/build/utils/run.php".
Nominatim and Postgres data are safely stored in a volume, so restarting the container does not erase data.


texts: Define polygon across the 180 meridian in DotSpatial

texts:
I am new to DotSpatial and I try to define a polygon that crosses the 180 meridian by following:

        FeatureSet fs = new FeatureSet(FeatureType.Polygon);
        Coordinate[] coord = new Coordinate[5];
        coord[0] = new Coordinate(30, 48);
        coord[1] = new Coordinate(120, 45);
        coord[2] = new Coordinate(-170, 64);
        coord[3] = new Coordinate(30, 70);
        coord[4] = new Coordinate(30, 48);
        Polygon pg = new Polygon(coord);
        pg.Normalize();
        fs.Features.Add(pg);
        fs.Projection = DotSpatial.Projections.KnownCoordinateSystems.Geographic.World.WGS1984;
        fs.SaveAs(HttpContext.Current.Server.MapPath("~/shp/test2.shp"), true);

However, when I check the shapefile, the polygon comes out to be

If I draw the same polygon in SQL Server with

DECLARE @Poly geography =
  geography::STPolyFromText('POLYGON((30 48, 120 45, -170 64, 30 70, 30 48))', 4326)
SELECT @Poly

I can see the polygon shown as

Can anyone help me how to solve this issue in DotSpatial?


texts: How to reproject a raster based on another raster when the datum is not said and get the same extent?

texts:
I am working with two rasters. I am trying to reproject r1 the same way as the raster r2. Rasters are available here

r1

class       : RasterLayer
dimensions  : 30286, 58025, 1757345150  (nrow, ncol, ncell)
resolution  : 99.99876, 99.99876  (x, y)
extent      : -2893269, 2909160, 224962.6, 3253525  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
values      : 0, 7 (min, max)

r2

class       : RasterLayer
dimensions  : 585, 1386, 810810  (nrow, ncol, ncell)
resolution  : 0.04166599, 0.041666  (x, y)
extent      : 235.207, 292.9561, 25.04224, 49.41686  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +no_defs
values      : -7.904284, 12.21496  (min, max)

I precise that the resolution of r1 is 100 m and the resolution of r2 is supposed to be 4km. Moreover, the system of coordinates is different as one is in degree and another one in meter.
I would like to reproject the rasters in order they keep their own resolution and share the same extent.

I did not find a way to do it. When I upload my data in Arcgis, the rasters are aligned as I would like which is not the case with Qgis and R. The final goal would be to allow raster calculations by stacking the raster for example.

I do not want to reproject the raster from r2 to r1 as I am working with thousands of raster that have the same coordinates than r2.

I can work with Arcgis, Qgis, and R. However, I do not know how to code in Python.

Currently, by applying the method by Spacedman (answers below), I obtained. This is much better as my rasters are closers but there are still not at the same position when I plot them. Moreover, the resolution of r1 has changed. I am confused to what is the problem.


texts: Using both leaflet's wms and extend on same script?

texts:
I have a main application which switches between multiple satellite options. Now the main application initially has all its options initialized using .wms like the code below:

              var satelitteLayer1 = L.tileLayer.wms(url,{
                   service: 'WMS'
                   //The rest of the header parameters
              });

now my challenge is there is another satellite API I would like to integrate but doesn't use .wms like above but .extend... When I add it to the same script I get an error message saying Map container already initialized, when I try adding it to the same as above I get an error message saying error authenticating although I did declare the key, id and authorization parameters. One of the tech guys from the company sent an example code which is instantiated like below:

             L.TileLayer.SecureLayer = L.TileLayer.extend({
                 createTile: function(coords,done){
                    //Logic of the function
                 }
             })

I am new in terms of using leaflet and would like to know if there a way to use both wms and extend on the same script since I see the approach is different? Or if anyone can actually explain the differences between .wms and .extend in layman's terms so that I can know which route to take in solving this...


texts: How to get something like psql within QGIS?

texts:
Does QGIS have something where I can type a PostGIS query, hit execute and it adds the results as a Layer in QGIS?

I realise I can create a view or table then open it in QGIS, but that becomes painstaking when you're trying out queries like you would with on the command line with psql.


texts: SelectFeature fails on a small polygon on top of a larger polygon

texts:
I'm using new OpenLayers.Control.DrawFeature to draw polygons and OpenLayers.Control.SelectFeature to select a specific polygon. I'm drawing one polgon and then a smaller polygon on top of the first one.

PolyA is a large polygon(50 x 50)
PolyB is a smaller polgon (10 x 10) drawn within (on top of) PolyA

Initially I can select and unselect PolyB correctly. However, once I have selected and unselected PolyA, PolyB is now hidden behind it and can no longer be selected.

Why is that and what can I do to resolve it?

select code (straightforward)

select  = new OpenLayers.Control.SelectFeature(layer,
     {
          highlightOnly: true,
          clickout: true, toggle: true,
          onSelect: function(feature) {

          //do stuff
     }
 });


texts:
I'm trying to download and process MODIS imagery using pymodis.  I have Python 2.7.8 (Spyder), numpy 1.9.1 and gdal 1.11.1.  I've also installed WXPython.

I originally tried to use the modis_download.py script, but continued to get

  SyntaxError: Invalid Syntax despite using the formatting recommended
  in the pymodis documentation:

modis_download.py -r -p MOD13Q1.005 -t h18v03 -f 2008-01-01 -e 2008-01-31 \modis_test

I'm now trying to call the pymodis.downmodis.downModis module to download MODIS data from the LP DAAC HTTP server.  This seems to be reaching the survey, but returns blank .txt and .log files:

pymodis.downmodis.downModis('E:\modis_test', url = 'http://e4ftI01.cr.usgs.gov', tiles = 'h26v08, h25v08', path = 'MOLT', product = 'MOD13Q1.005', today = '2014-01-01', enddate = '2014-02-02', debug = True)

Anyone familiar with pymodis that can lend a hand?


texts:
I have 2 issues:

With Qgis2web I can't export a heavy vector layer either as .shp or WMS.
I also can't make it display, the list of different values of the same attribute in the pop-up menu.


texts: Apply same symbology to many raster layers in ArcMap?

texts: I have many TIFF files (100+) that need to have the same symbology using the unique value option. Is there a quick way in arc 10.3 to do this?


texts: Making get feature request to forget between clicks in OpenLayers 3

texts:
I currently have a map with projects on it that when you click on it, it returns the information of said project using the code below.

The issue I'm having is that every so often it will duplicate the responses it is returning. It's almost as though it's remembering that it's been previously clicked on and then adding another response on top of it. How can I make it stop doing this?

  function projectName(browserEvent) {
    var coordinate = browserEvent.coordinate;
    var pixel = map.getPixelFromCoordinate(coordinate);
    var el = document.getElementById('projectNumber');
    el.innerHTML = '';
    map.forEachFeatureAtPixel(pixel, function(feature) {
      el.innerHTML += feature.get('Proj_Number') + '<br>';
      return true;
    });
  }
  map.on('pointerdown', projectName);

Edit: Tried testing this without the + in the += feature.get but it only returned one result when there are multiple points in the same place.


texts: How can I select which interpolation technique will be used in a gdal command? (gdal_fillnodata.py)

texts:
I work on Landsat 7 slc-off images and I would like to use the gdal_fillnodata.py command to fill-in the missing pixels. The problem that I have is that the gdal_fillnodata.py uses the "inverse distance weighting" (IDW) interpolation technique and I want to use the kriging geostatistical technique.

Do you have any idea how can I do it?


texts: Export 3D Shapefile Polygon to DWG as Polygon

texts:
I am looking for a way to export my 3D polygon shapefile to a .DWG polygon.  I have used the Conversion Tools - Export to CAD tool in ArcToolbox, but it will only export the polygon to line features.  Is there a setting or option that I am missing that will export as a polygon?

Thanks...


texts: Calculating distance between points of same layer using QGIS?

texts:
I am using QGIS 2.13 and a layer with hundreds of points. I need to calculate the distance to the nearest point... in the same layer. I tried the distance matrix but I get 0m since the "closest" point is the origin point.

The idea is calculate the average distance, std deviation, etc. and create a graph/histogram with these data to see how the distribution looks like.

Does anybody know how to calculate the distance to the nearest point in the same layer avoiding the origin point?


texts: ArcGIS Server security settings via script

texts:
Can anyone point me to a method for setting a service, folder, or the site root security setting to public via script? It can be using one of the admin tools (e.g. server\tools\admin\manageservice.py), arcpy, powershell, batch file, etc.

I have Desktop, Engine, ArcObjects, and WPF runtime SDK installed on the server. This is on Windows o/s.

My motivation is to be able to script/automate the server install and setup process, including services. I want to keep it as simple as possible because simple is easy to implement and debug. So if there is a solution that doesn't use the server's REST/SOAP interfaces, I would prefer to try that first.

However, I will consider REST if that is the best solution.


texts: mapserver: accessing mapfile with localhost instead of system path

texts:
Don't know if GIS SE is the right place to ask this question. Otherwise, please feel free to suggest me where should I migrate it.

I need to access my mapfile with localhost to test my service via my intranet.
At present, if I call:

http://localhost/cgi-bin/mapserver/mapserv.exe?map=C:/mapserver/umbe/mapfile/test7.map&mode=map

it works, while if I do:

http://localhost/cgi-bin/mapserver/mapserv.exe?map=http://localhost/mapserver/umbe/mapfile/test7.map&mode=map

it throws:

msLoadMap(): Unable to access file. (http://localhost/mapserver/umbe/mapfile/test7.map)

I guess it's a problem of folder/file permission, but I don't know how to solve it. I am using XAMPP, Apache web server with a Windows 10 OS. MapServer is version 7.0.4.


texts: GeoServer layer preview shows only missing images for tiles

texts:
Just getting started with GeoServer and already running into problems. When I click on Layer Preview and view any of the preloaded data in OpenLayers (or any other format) it is just blank and has the missing image icon where the data should be. One interesting thing is that I could click randomly in the blank map and see attributes of selected polygons appear so I think the data is partially loaded.
Example screenshot.

I've installed it on multiple machines and I am getting the same exact problem. I have connected my own data through PostGIS but I got the same result. That is why I am assuming it is user error. Anyone have any ideas?


texts: Why does ListLayoutElements give IndexError: list index out of range?

texts:
When trying to create a legend on ArcMap use the following Python/ArcPy code I received an error message:

import arcpy
mxd = arcpy.mapping.MapDocument(r"E:\cities\MyProject.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
lyr1 = arcpy.mapping.Layer(r"E:\cities\cities.shp")
lyr2 = arcpy.mapping.Layer(r"E:\cities\region.shp")
lyr3 = arcpy.mapping.Layer(r"E:\cities\lakes.shp")
lyr4 = arcpy.mapping.Layer(r"E:\cities\rivers.shp")
legend = arcpy.mapping.ListLayoutElements(mxd,"LEGEND_ELEMENT")[0]
legend.autoAdd = True
arcpy.mapping.AddLayer(df, lyr1, "BOTTOM")
arcpy.mapping.AddLayer(df, lyr2, "BOTTOM")
arcpy.mapping.AddLayer(df, lyr3, "BOTTOM")
arcpy.mapping.AddLayer(df, lyr4, "BOTTOM")
legend.adjustColumnCount(2)
mxd.save()

and I get this message :

Traceback (most recent call last):
  File "E:\cities\scripts\Lgend.py", line 8, in <module>
    legend = arcpy.mapping.ListLayoutElements(mxd,"LEGEND_ELEMENT")[0]
IndexError: list index out of range

What does this indicate?


texts: Can ArcGIS Desktop read any open source layer definitions?

texts:
I am looking for a non-proprietary way to store layer symbology for both vector and raster layers and then load those layers in ArcGIS Desktop 10.x.

I have used layer files for this purpose for many years. At runtime I load the layer file, retrieve the renderer and apply it to the layer I am loading into the current map document. Something like this:

IGxLayer pGXLayer = new GxLayer();
IGxFile pGXFile = (IGxFile)pGXLayer;
pGXFile.Path = symbology.FullName;

However layer files are proprietary and annoyingly version specific.

Is anyone aware of an open source symbology definition that ArcGIS Desktop can use? Something like QML? Is it true that ArcGIS Desktop cannot use QML files?


texts: How to Separate Bands Through ModelBuilder

texts:
Q1: I have a .tiff file containing 365 bands. Each band represents a daily variable. I want to separate each band and save as a single .tiff file. In the end, I want to have 365 individual .tiff files.

Q2: And since these are daily data, is there are a way to make them monthly if you can share your model.

How do I do this?

I am using ArcGIS 10.0.


texts: Attachment access using ESRI's dynamic services

texts:
I need to be able to display photo attachments in a tabbed popup window that ideally only contains a dynamic map service. I have a working demo running that uses a dynamic map to display the data and several feature services to gain access to the photo attachments and it works but I'd like to know if this is the only way to do it or if there's a simpler way...ie not having to use a Feature Service? The demo map only has two simple polygon's so it isn't that bad but my production map would have over 20 layers to loop though and there's a LOT of complex geometry so pushing that all out as a FS has proven to be much slower than a dynamic service.

The demo map I have running is available here: http://maps.ducks.ca/pmv

If do I have to utilize the FS method is there a better way than I've done here or is there someway to use the "showAttachments: True" option in the Popup properties?


texts:
How to add 3D features in Online ArcGIS (Story Map) - I created the 3D layers / features using ArcMap Desktop and ArcScene

I'm a student and a newbie looking for detailed information.


texts:
I think you know that the Mapzen service has closed. Unfortunately I'm on several projects based on R and Mapzen APIs, so I was looking for alternative services for which a package or a series of R functions were available (especially for geocoding).
Some geocoding packages I already considered:

opencage: I don't like the pricing and the provided features
tmaptools: there is a nice function, geocode_OSM, but since it uses OSM Nominatim servers after a certain amount of requests (i think ~2000) it is rate limited.
ggmap: I wouldn't use it since it uses Google services.

I think the best offer it's the Mapbox one (a lot of free requests and pay as you use) but i didn't find a r package for it yet.

Which other geocoding packages I should consider? (apart self hosting)


texts:
In QGIS 2.18, new vector layers resulting from QGIS tools have a default name related to the corresponding tool. For example, when I perform the "Lines to polygons" operation (under QGIS Vector menu), the resulting layer is always named Polygons from lines in the Layers Panel -- even if I ask for the shapefile to be saved under a specific name. When I compute a buffer, the new layer is called Buffer, and so on.

If I recall correctly, this was not the case in QGIS 2.14 (and previous versions): the default name of a layer was the name of its shapefile.

I find the new behaviour annoying since it leads to confusion between my layers; every time I use a QGIS tool, I have to uncheck the "Open output file after running algorithm" option, then manually open the newly created layer in order to display its "real" name. I can also change its name in the layer properties, but this is still irritating.

Is there a way to go back to the traditional behaviour? E.g. display by default the shapefile name (when specified) in the Layers Panel instead of some generic name?