texts: Can't set OGR layer geometry type in Python

texts:
I'm trying to create a layer of geometry type polygon, with a feature of geometry type polygon, and a geometry of type polygon. However, no matter how I attempt to set the geom type of the feature, geometry, or layer, the layer is always set to multipolygon.

Here is some example code to show what I mean:

# Create test polygon
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(1179091.1646903288, 712782.8838459781)
ring.AddPoint(1161053.0218226474, 667456.2684348812)
ring.AddPoint(1214704.933941905, 641092.8288590391)
ring.AddPoint(1228580.428455506, 682719.3123998424)
ring.AddPoint(1218405.0658121984, 721108.1805541387)
ring.AddPoint(1179091.1646903288, 712782.8838459781)
geom = ogr.Geometry(ogr.wkbPolygon)
geom.AddGeometry(ring)

print(geom.GetGeometryName()) # returns POLYGON

# Open the FGDB
driver = ogr.GetDriverByName("FileGDB")
path = "/path/to/polygon.gdb"
if not os.path.exists(path):
    data_source = driver.CreateDataSource(path)
else:
    data_source = driver.Open(path, 1)

# Define the spatial reference to use
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# Create the layer
layer_name = "polygon".encode("utf-8")
layer = data_source.CreateLayer(layer_name, srs, ogr.wkbPolygon)

print(ogr.GeometryTypeToName(layer.GetGeomType())) # returns Multi Polygon

# Create the feature
feat_def = ogr.FeatureDefn()
feat_def.SetGeomType(ogr.wkbPolygon)
feature = ogr.Feature(feat_def)

print(ogr.GeometryTypeToName(feature.GetDefnRef().GetGeomType())) # returns Polygon

# Set the geometry
feature.SetGeometry(geom)

print(feature.GetGeometryRef().GetGeometryName()) # returns POLYGON

# Create the feature in the layer
layer.CreateFeature(feature)

# Dereference the feature
feature = None

# Save and close the data source
data_source = None

The first print, I show that the geometry I created is of type
polygon
The second print, I show that the layer's geometry type is of
multipolygon, despite me initializing the layer with type polygon.
The third print shows that I did set the feature's geometry type to
polygon when explicitly creating a FeatureDefn
The fourth print shows that after setting the feature's geometry, it
is still type polygon

Here's the output of ogrinfo:

$ ogrinfo polygon.gdb/ polygon
INFO: Open of `polygon.gdb/'
      using driver `OpenFileGDB' successful.

Layer name: polygon
Geometry: Multi Polygon
Feature Count: 1
Extent: (1161053.021823, 641092.828859) - (1228580.428456, 721108.180554)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
FID Column = OBJECTID
Geometry Column = SHAPE
OGRFeature(polygon):1
  MULTIPOLYGON Z (((1179091.16469033 712782.883845978 0,1218405.0658122 721108.180554139 0,1228580.42845551 682719.312399842 0,1214704.93394191 641092.828859039 0,1161053.02182265 667456.268434881 0,1179091.16469033 712782.883845978 0)))

I've looked everywhere in the (somewhat useless) OGR Python bindings API and have found nothing useful. So, how do I actually set the layer to be of geometry type polygon?

This is with Python 2.7.15 and GDAL 2.3.2 on CentOS 7.5.


texts:
I would like to start using pure python to handle geospatial data since POSTGIS is not so optimized for rasters.

I starting using rasterio and with this code :

with rasterio.open(clc_r) as src:
    print(src.width, src.height)
    print(src.crs)
    print(src.transform)
    print(src.count)
    print(src.indexes)
    print(src.profile)
    array = src.read()
    array.shape

I got this output:

63976 45242
+ellps=GRS80 +lat_0=52 +lon_0=10 +no_defs +proj=laea +units=m +x_0=4321000 +y_0=3210000
| 100.00, 0.00, 918900.00|
| 0.00,-100.00, 5440600.00|
| 0.00, 0.00, 1.00|
1
(1,)
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': 65536.0, 'width': 63976, 'height': 45242, 'count': 1, 'crs': CRS({'proj': 'laea', 'lat_0': 52, 'lon_0': 10, 'x_0': 4321000, 'y_0': 3210000, 'ellps': 'GRS80', 'units': 'm', 'no_defs': True}), 'transform': Affine(100.0, 0.0, 918900.0,
       0.0, -100.0, 5440600.000000001), 'tiled': False, 'compress': 'deflate', 'interleave': 'band'}
Traceback (most recent call last):
  File "raster_calculation.py", line 23, in <module>
    array = src.read()
  File "rasterio/_io.pyx", line 316, in rasterio._io.DatasetReaderBase.read
MemoryError

I am getting crazy trying to understand why with:

RAM 7,8 GiB
 os 64-bit

I am not able to read a raster of 14 MB... What am I doing wrong?
If even with such a small raster I got an issue like that, the processing will be even more tricky. (and with larger files even worst..)

I am running on Ubuntu 18.04 using python 3 64 bit..

UPDATE

The file I was trying to read had a lot of NoData compressed with the LZW compression (i dunno much about this: just figured it out ...).
I didn't notice this on QGIS.. So it's the extension of the entire Europe!
I clipped it better and it is ok now!


texts: Getting event object from Draw Toolbar onDrawEnd using ArcGIS API for JavaScript?

texts:
I am using the Draw Toolbar in the Esri Javascript API v3.4 and connecting to the onDrawEnd event using aspect.after(drawingToolbar, "onDrawEnd", drawEndHandleFunction, true);

I am using the Draw Toolbar with the EXTENT tool to allow my user to select features by drawing a box around them.  So far it is working fine.  However, I want the user to be able to add and subtract features from the selection using the CTRL and SHIFT keys while drawing.

I was able to do this in Chrome using the event.ctrlKey and event.shiftKey properties of the event global variable.  Firefox does not appear to have an event global variable like this and so my function doesn't work.  Firefox apparently passes an event as a parameter in the function, but that is not available because I used dojo/aspect.

Is there a better way to handle this?  Otherwise is there a way to get a map mouse event that triggers the onDrawEnd event in the Draw Toolbar?


texts: Auto update map service when feature class is updated

texts:
I wrote a python script to create an MXD from template and auto publish it.

The following is the excerpt of my code:

#publish MXD method
def PublishMXDToArcGISServer(id, in_arcgisAgs, in_publishWS):
  sddraft = "{}\\SDDRAFT\\{}.sddraft".format(in_publishWS, id);
  service_name = id.replace('-', '_').replace('.','')

  # add the layer file to your map
  mxd_location = "{}\\MXD\\{}.mxd".format(in_publishWS, id)
  mxd = arcpy.mapping.MapDocument(mxd_location)

  #Create SDDRAFT file
  analysis = arcpy.mapping.CreateMapSDDraft(mxd, sddraft, service_name, 'ARCGIS_SERVER', in_arcgisAgs, False, "Folder");

  sd = "{}\\SDDRAFT\\{}.sd".format(in_publishWS, id);

  arcpy.StageService_server(sddraft, sd)

  arcpy.UploadServiceDefinition_server(sd, in_arcgisAgs)

Above code can successfully published the map service, but when I added a new field or alter a field in the feature class, the updated field is not reflected on the Map Service.

But when I used an ArcGIS Desktop to publish a map service, any new field or altered field in the under lying feature class, is reflected on the map service.

Is there any thing i need to specify when publishing using ArcPy in order for the updated field to be reflected on the Map Service?

Updating Web Map Service (WMS) by default when data amended in Oracle source tables? mentioned a similar issue but I do not understand the answer given.


texts:
How big is a MODIS 250m pixel in reality?

This appears to be a question that contains its own answer, but I have several reasons to doubt that the pixel size is actually 250m x 250m:

I downloaded MODIS 250m NDVI data from the USGS site, transformed it from Sin to Geographic (SRID 4326) using the MODIS Reprojection Tool (which has since been replaced, but at the time was the official tool).
Using gdalinfo on the resulting file gives the pixel size as (0.002884053983564 X 0.002884053983564).  I assume that this measurement is in decimal degrees.
This online calculator gives the length of one degree of latitude at the latitude of my study area as 111,092.7 m, and of longitude as 81,540.9 m.  That would make the pixel dimensions 320m in the N-S direction and 235m in the E-W direction.
I've overlaid the NDVI image onto other data, and the two correspond; for example you can clearly see greener vegetation along rivers.  This wouldn't happen if they were misaligned or mis-projected.
I can also measure the size of a pixel using QGIS's measurement tool, and the length comes out to about 320 x 235 m, i.e., it agrees with the file metadata.

I have read the documentation and found a few suspicious references in published papers to things like "nominal" pixel size of 250 m, or pixel size of 250 m "at nadir", and this paper makes it look to me as if the physical camera/mirror apparatus would result in pixels that aren't 250m x 250m square, but the authors don't describe in detail the projection and/or transformation used to create the data products so I'm not sure.


texts: v.what.rast cannot find points map layer

texts:
I would like to apply v.what.rast with a particular points map. Trying the straightforward way, the module complains there is no layer 1 (the default) in the points map:

> v.what.rast map=debug_interpol_points_corr raster=weights column='value'
ERROR: Database connection not defined for layer 1

First the sanity check, there is indeed a layer 1 in this map:

> v.category input=debug_interpol_points_corr option=report
Layer: 1
type       count        min        max
point         30          1         30
line           0          0          0
boundary       0          0          0
centroid       0          0          0
area           0          0          0
face           0          0          0
kernel         0          0          0
all           30          1         30

Trying a different approach, with the -1 layer number to force the module to use all the layers in the map:

> v.what.rast map=debug_interpol_points_corr raster=weights column='value' layer=-1
ERROR: Database connection not defined for layer -1

What is this module exactly looking for? Or is this a bug?

Update: This is the output from v.db.connect:

> v.db.connect -p map=debug_interpol_points_corr
Vector map <debug_interpol_points_corr> is connected by:
layer <2/debug_interpol_area_points_2> table <debug_interpol_points_corr> in database </home/desouslu/EAWAG/Projects/RasterStitching/grass/Luzern/Merging03/sqlite/sqlite.db> through driver <sqlite> with key <cat>

Is there anything wrong with it?


texts: Set projection for all features generated in PYQGIS script

texts:
I have a script that loops through input raster's creates new raster files and converts these to vector performing a number of additional processes. Ultimately a vector is created that is a modified version of the input raster.

My question is how can I ensure that each new raster and vector created, whether they are a temp file or the final output, maintain the same projection (ESPG 27700)? Is there a command I can add in at the start of each loop or can I set the output projection in each individual process e.g. raster calculator or polygonize. An example line of the script is provide:

out1=processing.runalg('gdalogr:rastercalculator', raster,'1',None,'1',None,'1',None,'1',None,'1',None,'1','A>0',None,5,None,result_shp + "/" + result_name)

Additionally, is it possible to set a snap raster in a similar way and is this something I should consider doing?


texts: Problems rendering tilestache served mbtiles in QGIS

texts:
I came across Geoffery Hing's walk through on using Tilestache to serve up mbtiles to QGIS. I followed the configuration and I was able to view tiles in the browser without any problems. However, when I try to load the XML in QGIS, all I get is an image that looks roughly like scrambled cable from the days of yore:

My tilestache.xml is as follows:

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>http://127.0.0.1:8080/${z}/${x}/${y}.png</ServerUrl>
    </Service>
    <DataWindow>
        <UpperLeftX>-73.009</UpperLeftX>
        <UpperLeftY>41.3578</UpperLeftY>
        <LowerRightX>-72.8394</LowerRightX>
        <LowerRightY>41.2355</LowerRightY>
        <TileLevel>13</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:900913</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <Cache />
</GDAL_WMS>

I have a hunch that this is a problem with tilestache and perhaps its dependencies. GDAL is version 1.9, which should be sufficient. I have no problem loading tiles from openstreet map using this XML.

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>http://tile.openstreetmap.org/${z}/${x}/${y}.png</ServerUrl>
    </Service>
    <DataWindow>
        <UpperLeftX>-20037508.34</UpperLeftX>
        <UpperLeftY>20037508.34</UpperLeftY>
        <LowerRightX>20037508.34</LowerRightX>
        <LowerRightY>-20037508.34</LowerRightY>
        <TileLevel>18</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:900913</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <Cache />
</GDAL_WMS>

When I run gdalinfo on the mbtiles, I get an error saying that they are not a recognized file format. I suspect this is at the root of my problem, but I am at a loss as to how to fix this.


texts: ArcGIS 10.2 Calculate xmin, xmax, ymin, ymax polygons

texts: In ArcGIS 10.2, I would like to know how can I calculate xmin, xmax, ymin, ymax coordinates in a polygon features?


texts: How does QGIS open so large raster datasets (about 40GB)?

texts: I have problem with GDAL library when open a big GeoTiff file, with size about 32000x32000. I cannot use ReadAsArray function because maximum size of numpy array in python. But I'm wondering why QGIS can open that file easily. What is the technique behind?


texts: MarkerSymbolizer of MultiLine not orientated in line direction (Mapnik, Mapbox)

texts:
What I want:

drawing orientated arrows on one-directional streets

What I tried so far:

MarkerSymbolizer: However, the arrows are all looking towards the "east" (right of the screen) instead of the direction of the MultiLineGeometry

Code:

<Rule>
  <LineSymbolizer />
  <MarkersSymbolizer marker-placement="line" file="arrow.png" width="1" />
</Rule>

The data which is rendered here is a MultiLineString.

UPDATE:

The only thing which I am missing right now in my solution posted below (using LinePatternSymbolizer) is that I cannot define spacing and max-error etc like here. Does anyone know if this is possible?


texts:
I have an aerial photo which consists of 2 files (.jp2, j2w) and i want to perform a segmentation.

To be exact, i can't upload the 2 files with r.in.gdal because it says that the data file cannot be supported but in order to do segmentation (i.segment) it requires to do a group  with these files, with i.group.

The question is what to do with the files in order to do group and segmentation.


texts:
I am developing a data import script tool for ArcGIS that downloads from an SQL Server into a pandas DataFrame. This python script worked perfectly as I built it in Spyder however now that I am setting it up in a Toolbox, ArcGIS cannot import pandas.

I have tried Installing pandas with Esri's python? on my own however what I am looking to do is automate the process so as soon as the tool is run the script checks for pandas and installs as necessary on any user's machine. I already fixed this similar problem with pypyodbc by storing the necessary setup scripts in the public drive folder where the tool script is located:

import os
import imp
path = os.path.realpath(__file__)
path = path[:-16]
os.chdir(path)
try:
    imp.find_module('pypyodbc')
    found = True
except ImportError:
    found = False
if found == False:
        import setup
        import pypyodbc

I can't seem to apply this method to pandas as there isn't an install/setup script I can save in my public folder. Is there any way to automate the install of pandas from a Script Tool in ArcGIS?


texts: Displaying coordinates and inputs as LatLon or LonLat?

texts:
I'm trying to get a sense if this is an issue for others or every input/output should be labeled so user is not confused and just go with it?

I think almost everyone pronounces it as "LatLon".

Who started it?

Is it because it's in alphabetical order compared to "LonLat"?

Mapping Lat and Lon to Cartesian plane Lon is "x" and Lat is "y" so since we say "(x,y)" it should be said as "LonLat".  And now for display of information.

Should the status bar on a mapping application display La,Lo or Lo,Lat?

Should it just be labeled as one way and let user deal with it?

And same with input, what's the right way to order the fields?

KML's format is Lon,Lat,Altitude.  While other apps is Lat,Lon and so have to be very vigilant about when converting formats.

Is there a standard?


texts: OpenLayers 4.x DescribeFeatureType to populate a filter menu

texts:
With OpenLayers 4, I'm trying to get field names from a WFS request (https://localhost:8080/myws?SERVICE=WFS&REQUEST=DescribeFeatureType&VERSION=2.0.0&TYPENAME=ms:layername&outputFormat=gml/2.1.2) and parse it into JSON or something else.

I found a solution, but it is for OpenLayers 2 (example)

Does somebody know if

OpenLayers.Format.WFSDescribeFeatureType();

will be deployed into OpenLayers 4?

Any other way possible? My goal is to "harvest" my WFS DescribeFeatureType to populate a filtering menu.


texts: Error: Adding WCS ( Web coverage Service ) into ArcMap 10.4.1

texts:
I am trying to publish some raster data sets as WCS using ArcGIS 10.4.1 (Desktop ArcMap and Server Manager), but have had little success. I've attempted to publish a data set stored as an ESRI Grid, a Geo-tiff, and a Mosaic data set. I am able to connect to the WCS image server, but when I try to add the image I get the following message:

"Could not add the specified data object to the map. Failed to open raster dataset."

Here are detailed steps:

First I created the shared folder to register my data source with ArcGIS server. Next in ArcMap Catalog I navigated to to the registered folder and selected the required raster data set ( Geo-tiff both geographic and sinusoidal ) also created a Geo-database ".gdb" for a list of mosaic data set.

Next using ArcMap publishing Wizard, I have published the Image services as WCS and WMS both. When publishing my dataset, I select the WCS and WMS capabilities along with the Image Service and verified on the Server Manager that the services I created are available on server along with under WMS and WCS tabs.

Next I needed to check if I can get these WMS and WCS links work in ArcMap. I am able to add the image as a Image Service when I connect to the ArcGIS Server just fine. I'm also able to add the published data as a WMS ( can see the image in Data view ). The WCS is the only service that fails. Does anyone know what would cause my WCS to fail, but allow my Image Service and WMS work.

There was a bug report of dragging the WCS layers in Data view of ArcMap might be issue so I tried the add layer button to select the WCS layers from WCS server connection but it still failed.

Below link, it seems like when we add WCS service in ArcMap it atleast gives a black layer which can be contrast stretched to get the right image. But in my case I am not getting any layer just the error mentioned above.

http://docs.gina.alaska.edu/ndvi/how_to.html

Also there seems to be not issue in my Geo-tiff's as when I add these files directly to ArcMap it will show up layers / images also using WMS it works.

Last, I also tried the same WCS service GET Coverage operation in browser and I got a Tiff Image downloaded. ( So same WCS service working in browser but not in ArcMap )

https://localhost:6443/arcgis/services/imageservicename/ImageServer/WCSServer?SERVICE=WCS&VERSION=1.1.1&REQUEST=GetCoverage&IDENTIFIER=1&FORMAT=image/GeoTIFF&BOUNDINGBOX=-124.45790612602009,40.000531011912074,-91.378950375060455,49.999999995506855,urn:ogc:def:crs:EPSG::4326


texts:
There is a task before me of generating images as shown below. The main problem is the signs which must be placed outside the polygon (if they do not fit in it), and the lines should connect the sign with the point on the polygon.

Are there any solutions exist in free map renders?

Map example:


texts: What effect does standard deviation have on the display of DEM rasters?

texts:
I need to explain why ArcGIS offers the function "standard deviation" when symbolizing rasters with the Stretch method.

The Standard Deviation is a measure of how spread out numbers are and is calculated by the square root of the Variance. The variance in the Squared differences from the mean and the "n" in the ArcGIS Setting refers to the population.

So what is really happening? Is ArcGIS replacing each Height value with a StdDev value calculated by neighbouring "n" pixels? That would mean the colours represent the average "spread" of heights and not the actual height values themselves.

Am I on the right track here?


texts: Returned non-zero exit status 2 gdal_merge.py in QGIS plugin

texts:
I'm trying to merge several files within a QGIS plugin via gdal_merge.py. In Anaconda distribution exactly the same code works as it is like that. I suppose Python in QGIS can't find gdal_merge.py. But maybe I'm also wrong. How can I fix this issue?
wd = C:/Users/Jonas/Documents/Bachelorarbeit/MODIS/Data/oneday
# Define directories
Path(wd+"/mergedTiles").mkdir(parents=True, exist_ok=True)
out_dir = wd + '/mergedTiles/'

# Get tif files in wd
tile_list = glob.glob(wd+'/unmergedTiles' + '/*LST_*.tiff')

#get sublsists and zip
str_start = [s.split('\\')[1][:17] for s in tile_list]
str_end = [s.split('\\')[1].split('_')[1] for s in tile_list]
res = np.unique([i + j for i, j in zip(str_start, str_end)])

for label in res:
    files = [wd+'/unmergedTiles/' + s.split('\\')[1] for s in tile_list if s.endswith(label[17:]) and s.split('\\')[1][:17] == label[:17]]
    out_path = out_dir + label
    #print('Files To Merge')
    #print(files)
    cmd = 'python ' + "gdal_merge.py" +' -co COMPRESSED=YES'+ ' -o' + ' ' + out_path + ' ' + '-of gtiff ' + " ".join(files)
    #print(cmd)
    output = subprocess.check_output(cmd)


texts:
In arc gis for server 10.2 I have deleted some map services. But in the server logs it shows the following messages regularly even though the services are not in use anymore. I deleted the services through server manager

SEVERE

29 Apr 2015 13:25:08

Unable to process request. Service (service name would be here) not found

Rest

WARNING

29 Apr 2015 13:25:08

Error getting service.

Rest

Could anyone suggest why it may be logging these errors regularly?


texts: unselect all Features with jquery dialog close button openlayers 2.13

texts:
how can one unselect all Features with jquery dialog close button function in openlayers 2.13?

View Entire Project Here (some buildings are not yet clickable)

incorrect code:

OpenLayers.layer.selectedFeatures.unselectAll()

code in relation to function:

        var selectControlClicks = new OpenLayers.Control.SelectFeature(vectorLayer, {
        onSelect: function(feature) {
            $.get( "infoget.php", { name: feature.attributes.name } )
                .done(function( data ) {
                    var options = {
                        $FillMode: 1,
                        $DragOrientation: 1,
                        $ArrowNavigatorOptions: {
                            $Class: $JssorArrowNavigator$,
                            $ChanceToShow: 2,
                            $AutoCenter: 0,
                            $Steps: 1
                        }
                    };
                    dialog = $( data ).dialog({ minHeight: 390,maxHeight: 390,draggable: false,close: function( event, ui ) {OpenLayers.layer.selectedFeatures.unselectAll()}} ).tabs();
                    var jssor_slider1 = new $JssorSlider$("slider1_container", options);
                });
        },
        onUnselect: function(feature) {
            dialog.dialog("destroy").remove();
        },
    });


texts: QGIS styling - displaying an outline around line features in QGIS

texts:
I am trying to symbolize a road network in QGIS.

Currently I am using a pen width of 1.86 milimeters to symbolize the roads, however I would also like to outline the roads with a secondary color.

This is easily done for polygons, however I cannot figure out if it is also possible with line features.

Another solution would be to use a buffer to generate polygons from the roads, however this is not optimal, so I would like to avoid this.

FINAL EDIT: I do not think this question exactly has been answered before, at least I couldn't find it. However, the answers in the linked questions could also be applied for my case.


texts: Why are values different between attribute table in QGIS and dbf?

texts:
I'm working on a shapefile on QGIS 2.8.6, on Debian Jessy. Its attribute table contains two columns (real 8,6) for longitude and latitude coordinates. When working from QGIS I'm able to write normal values, i.e :

longitude;latitude
1.153303;43.834765
1.136567;43.822574

However when I try to open the dbf file with LibreOffice Calc, strange things occur : the data type is set to "N,9,6" instead of 8,6 and about 80% of the values are rounded, i.e. :

longitude;latitude
1.000000;43.000000
0.000000;43.000000

This phenomenon is the same when I try with OpenOffice or Microsoft Excel (with one exception : with OpenOffice I was able to get the correct values but only once).

Same goes when I try to save the shapefile from QGIS in a CSV file ; values get "lost in translation".

This seems totally random to me, as I'm unable to explain this behaviour. I tried copying values from longitude and latitude columns in text columns in QGIS, and then opening the dbf : in this case the original longitude and latitude all appear with the correct value (however with real 9,6 and 7,6 this time!).

Does someone have an explanation ? Does this come from QGIS itself or is this a problem when reading the dbf file ?

Edit : I wonder if the problem could come from the decimal mark ? I used the comma which is the decimal mark used in my country but it seems that some softwares are waiting points. Could it be a reason ?

Edit2 : I found a way to export a CSV of the attribute table and getting the correct values, by using the MMQGIS plugin (as proposed in this answer), which allows you to explicitely state the decimal mark. It works where the native option in QGIS didn't.

This is not an explanation nor a solution to my initial problem though, but I guess it could be useful to some people.


texts: Vector feature "flashes" after geometry update in OL-Cesium

texts:
I would like to append new points to a vector polygon geometry in OpenLayers v6.5.0 and Ol-Cesium v2.12 (Cesium v1.77). My UpdatePolygon function will be called in every seconds or more frequently with new parameters. In 2D mode (if the cesViewer is disabled) it works perfectly. However, in 3D mode (cesViewer enabled) the polygon flashes after every update. It is annoying because the polygon will be updated in every second or more frequently and the view is following the polygon, thus it seem like flashing continuously. The "camera" height is very low (zoom level is the highest), and the view is "moving" continuously, like a car navigation in 3D view.
The user constantly sees the polygon. The feature is a single, simple polygon with only one, exterior ring.
This is the relevant part of my code:
function updatePolygon(x_prev,y_prev,x,y,polygonLeftSideDistance,polygonRightSideDistance){

heading = turf.bearing(prev_position, current_position);

//Calculate new coordinates
var currentGeom = source.getFeatureById(currentID).getGeometry().getCoordinates()[0];
currentGeom.pop();
currentGeom.push(turf.toMercator(turf.destination(current_position,polygonLeftSideDistance,heading-90)).geometry.coordinates);
currentGeom.unshift(turf.toMercator(turf.destination(current_position,polygonRightSideDistance,heading-90)).geometry.coordinates);
currentGeom.push(currentGeom[0]);

//Update the geometry
source.getFeatureById(currentID).getGeometry().setCoordinates([currentGeom]);
}

//Source
var source = new ol.source.Vector( {overlaps: false,strategy: ol.loadingstrategy.bbox});
var format = new ol.format.GeoJSON();
//Vector layer definition
var vectorLayer = new ol.layer.Vector({
        renderMode: 'image',
        source: source,
        disableHitDetection: true,
        style:  function(feature, resolution) {
    const name = feature.get('category');
    return name == "1" ? style1 : style2;
  }
 });

//Map Creation
        var map = new ol.Map({
        renderer: 'webgl',
            controls: controls,
            interactions: interactions,
            updateWhileInteracting: true,
            updateWhileAnimating: true,
            target: 'olCesMap',
            layers: [basemapTileLayer, vectorLayer],
            view: new ol.View({
                center: center,
                zoom: zoom
            })
        });
        //Cesium Viewer constructed based on Ol.Map
        var cesViewer = new olcs.OLCesium({
            map: map,
            sceneOptions: {
                mapProjection: new Cesium.WebMercatorProjection()
            }
        });
        //Cesium Viewer scene
        var scene = cesViewer.getCesiumScene();
        //Cesium Viewer enable
        cesViewer.setEnabled(false);

        //Cesium Viewer depthtest
        scene.globe.depthTestAgainstTerrain = false;

         function setCesViewer3D(){
        if(!cesViewer.getEnabled()){
            cesViewer.setEnabled(true);
            cesViewer.scene_.debugShowFramesPerSecond = true;
cesViewer.scene_.requestRenderMode = true;cesViewer.enableAutoRenderLoop();
        }

Is it a bug in OL-Cesium? Do you have any suggestions how to fix it?


texts: Using Extent and SpatialReference with ArcGIS API for JavaScript?

texts:
I am trying to display an ArcGIS Map on my webpage.

Whenever I try to use the 'extent' Attribute next to the 'center' and 'zoom' Attributes, I get the following error:

GET https://js.arcgis.com/4.8/esri/SpatialReference.js net::ERR_ABORTED 404

So apparently the SpatialReference.js can not be loaded, but is a part of the ArcGIS API.

Am I misunderstanding something?

I have my code from https://developers.arcgis.com/javascript/latest/api-reference/esri-views-MapView.html#using where the use of 'extent' is explained.

<code>
<script src="https://js.arcgis.com/4.8/"></script>
// ...
require([
            "esri/views/MapView",
            "esri/geometry/Extent",
            "esri/SpatialReference",
            "esri/WebMap",
            "dojo/domReady!"
        ], function (MapView, WebMap, Extent, SpatialReference) {
            var webmap = new WebMap({portalItem: { id: "f7f....."}});
            var view = new MapView({
                map: webmap,
                container: "mapdiv",
                center: [7.16809, 51.94349],
                zoom: 14
            });
            // Set the extent on the view
            view.extent = new Extent({
                xmin: -9177882,
                ymin: 4246761,
                xmax: -9176720,
                ymax: 4247967,
                spatialReference: new SpatialReference({ wkid: 102667 })
            });
        });
</code>


texts: Spatial ETL Tools won't run "AttributeError: 'module' object has no attribute, etc."

texts:
I have a python script that runs 5 Spatial ETL tools that I wish to run as a scheduled task.

Each of the tools will run successfully from ArcCatalog and also from the ArcCatalog Command line.

If I try to run the tools from the DOS window, as would be done via python script for a scheduled task, arcpy cannot find the tools, although it does load the toolbox.

In the attached screenshot:

I open a cmd.exe window
start python
import arcpy
Import my toolbox, which is verified as loaded on the next line
Use arcpy.ListTools to see what tools are loaded. the result is "[]", meaning no tools are recognized from this tool box
I use the same ListTools to show the Spatial Analyst tools that were loaded by import arcpy
I attempt to execute one of the tools from my toolbox and forget a " (my mistake!)
I attempt to execute one of the tools from my toolbox and receive the message "AttributeError: 'module' object has no attribute 'CopyFeatureClass_newgpetl'".

What is preventing my tools from being recognized from the DOS CMD.EXE window while ArcCatalog and ArcCatalog Command Line can run them?


texts: I've made a custom plugin and use QT Designer. When the form appears, if I click into the cavenas then the form disappears. I would like to force the forme to stay a the first plan. And also, keep the "possibility" to click into the cavenas if I need (without hiding the plugin's form). How can I manage that?


texts:
Is there any possibilty to add a little WindowsForms with a textbox to an AddIn-Toolbar? In ArcGIS 9.3 I could easily add it to the Toolbar.cs like all the buttons ect with AddItem. But how can I do this in ArcGIS10?

To say it in more details:
What I want to do is looping throught features. Of course I have a NextFeature- and a PreviousFeature-Button. But I also want to give the possibility to go to a certain feature by entering the OID in a textbox.

I'm thankful for any suggestion.


texts: Access features with click events in SceneView ArcGIS Javascript API?

texts:
I am trying to access features of a featurelayer in a SceneView on click event. I have seen this example from Esri's documentation for ArcGIS API for JavaScript but this example is for a MapView. When I change MapView to SceneView, the graphic property of the object that is returned by view.hitTest is always null, even if I click on a feature.

I simply just replaced every MapView in the code to SceneView. The map loads fine, click event works (i.e. it returns x,y coordinates of the click).


texts: Using dot notation vs underscores when calling ArcPy functions?

texts:
I'm working on a Python script that is heavily utilizing Network Analyst tools. I understand that dot notation is used to call methods from a class, but don't really get why/how ArcPy functions appear to be callable with different syntaxes. Example:

     import arcpy
     from arcpy import na

     # what is the difference between this:
     arcpy.na.CopyTraversedSourceFeatures(routeLyr, gdb, edges, junctions, turns)
     #and this:
     arcpy.CopyTraversedSourceFeatures_na(routeLyr, gdb, edges, junctions, turns)

Both syntaxes appear to behave the same when running from the Python window in an ArcMap session, but will these behave differently when running from a different environment, such as a script tool or a batch file run in Task Scheduler?


texts:
I have created the web application using the default setting of QGIS2Web for OpenLayers but now want to customize the tools, e.g. increasing Search and Zoom button size and appearance, adding Pan and Fit to display option.
Where can I add these changes in code and in which file?


texts: Exporting dataframe - used reduceRegion to sample an image in Google Earth Engine

texts:
I am trying to sample the Global Forest Change dataset using randomly drawn points and the reduceRegion function. The way I construct this is:

var sampledPoints = gfc2014_c.reduceRegions({
  collection:points,
  reducer: ee.Reducer.mean()
});

Where gfc2014_c is the GFC datset cropped to sub-Saharan Africa and points are my randomly drawn points. The entire code is here.

I was aiming for a large sample of 1 million pixels. When I try to export in chunks (even chunks as small as 5,000) I get the following error:

Error: image.reduceRegions: Computed value is too large.

What am I doing wrong? I know the GFC dataset is really large, but it must be possible, right?


texts:
I have been trying to merge about 40 layers for a while now but I keep getting this error message:
{} field in layer {} has different data type than in other layers.

I have actually done in before and it worked and I have to redo it again but for some reason it is not working. Any help please on what the problem could be?


texts: PyQGIS 3 Single user shapefile editing for network drive

texts: I have a project saved on a shared network drive, with a group of people working on it. I need a PyQGIS solution to prevent multiple users from editing the same shapefile. That is, once one of the users has switched a layer into 'edit mode', other users should not be able to do the same until they have disabled edit mode and saved their changes, so as not to overwrite each others' changes. A simple script would suffice, I will then code it into a QGIS plugin.


texts:
I am brand new to QGIS.

I have downloaded the BLM's public land survey system from their website (https://catalog.data.gov/dataset/blm-national-public-land-survey-system-polygons), which has a file format of .lpk.

I cannot seem to open that file format in QGIS.

Anyone know how to get around this problem?


texts:
I have some data in a CSV file. One of the interested variable is test done (1=done, 2=not done). There is a variable of county name in which I would like to join with a county shapefile.

I only know how to join them using the county file as the base using ArcGIS, but the interested variable cannot be summarized correctly (as I would like each county to get a percent of 'done' from the interested variable).

How do I join it so I use the CSV as a base, such that I can get a percentage of 'done' later on?


texts:
I recently came across a loop_zonal_stats function by @ustroetz found here: Issue Trying to create Zonal Statistics using Gdal and Python.. I am using a (single) polygon shapefile and a raster file to compute the zonal statistics. The polygon shapefile (http://www.epa.gov/wed/pages/ecoregions/na_eco.htm#Downloads) is a list of North American ecoregions. I first opened the file in QGIS and saved it in the EPSG:4326 - WGS 84 coordinate reference system. The raster file (http://worldclim.org/current) contains precipitation data and is an ESRI grid file with 30 arc-second resolution.

One issue I noticed with the code is that NoData values are included in the calculations (a "bug" that I noticed in QGIS as well). In order to fix this, I added the following lines of code to the zonal_stats function:

nodata = banddataraster.GetNoDataValue()
if nodata is not None:
    zoneraster = np.ma.masked_equal(zoneraster, nodata)

After making this change, I am successfully able to exclude the NoData values from the calculation. However, I noticed something strange that I don't understand.
Many of the polygons in my shapefile have the same attribute (and none of them overlap). I used QGIS to convert the shapefile of polygons to a shapefile of multipolygons (QGIS SinglePart to MultiPart) and then applied the loop_zonal_stats function to compute the sum and count for this multipolygon shapefile. However, the results from the original polygon shapefile are inconsistent with those from the multipolygon shapefile. For example, for a particular region with "NA_L3CODE" = 10.1.2, I get a count of 352 + 6873 + 317400 = 324625 for the original polygon shapefile and a count of 369288 for the multipolygon shapefile. Might anyone know what is going on?


texts: QGis-Server 2.14 raises error reading project file on Ubuntu 16.04

texts:
I installed QGis-Server 2.14 on Ubuntu 16.04 following the instructions on this blog. The fgci on its own is responsive and correctly replying to a GetCapabilities request like:

http://localhost/cgi-bin/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

However, if I point the browser to a QGis project file, like:

http://localhost/cgi-bin/project/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

The response is an error:

<ServiceExceptionReport version="1.3.0">
  <ServiceException code="WMS configuration error">
    There was an error reading the project file or the SLD configuration
  </ServiceException>
</ServiceExceptionReport>

Using directly the command line I get a different error message:

$ ./qgis_mapserv.fcgi "QUERY_STRING=SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities"
QFSFileEngine::open: No file name specified
Warning 1: Unable to find driver ECW to unload from GDAL_SKIP environment variable.
Warning 1: Unable to find driver ECW to unload from GDAL_SKIP environment variable.
Warning 1: Unable to find driver JP2ECW to unload from GDAL_SKIP environment variable.
Warning 1: Unable to find driver ECW to unload from GDAL_SKIP environment variable.
Warning 1: Unable to find driver JP2ECW to unload from GDAL_SKIP environment variable.
Content-Length: 206
Content-Type: text/xml; charset=utf-8

<ServiceExceptionReport version="1.3.0" xmlns="http://www.opengis.net/ogc">
 <ServiceException code="Service configuration error">Service unknown or unsupported</ServiceException>
</ServiceExceptionReport>

The Apache logs are clean. The error log is empty and the access log is reporting nothing unusual.

Are there any other logs I can look into? Any hints on what may be causing this error(s)?

Update I: I was able to set up logging following the suggestions by elpaso66. As he expected, QGis-Server is returning a file not found error, but the file actually exists:

$ tail /tmp/qgis-000.log
[8066][15:50:44] inserting pair MAP // /usr/lib/cgi-bin/project/project.qgs into the parameter map
[8066][15:50:44] MAP:/usr/lib/cgi-bin/project/project.qgs
[8066][15:50:44] REQUEST:GetCapabilities
[8066][15:50:44] SERVICE:WMS
[8066][15:50:44] VERSION:1.3.0
[8066][15:50:44] Error, configuration file '/usr/lib/cgi-bin/project/project.qgs' does not exist
[8066][15:50:44] Checking byte array is ok to set...
[8066][15:50:44] Byte array looks good, setting response...
[8066][15:50:44] Sending HTTP response
[8066][15:50:44] Request finished in 0 ms

$ ls -la /usr/lib/cgi-bin/project/project.qgs
lrwxrwxrwx 1 root root 49 Sep 29 14:41 /usr/lib/cgi-bin/project/project.qgs -> /home/desouslu/Projects/project/project.qgs

$ ls -la /home/desouslu/Projects/project/project.qgs
-rw-rw-r-- 1 desouslu desouslu 156185 Aug 16 15:40 /home/desouslu/Projects/project/project.qgs

Update II: This ended up being an issue with permissions. On Ubuntu, folders created in the home folder by default do not have execution permissions for all. A chmod a+x is required on the folder structure leading to the .qgs file, otherwise Apache fails to access it.


texts: Efficiently query two properties QGIS spatialite

texts:
I have a simple query that involves two steps:

non-spatial criteria (should be super-fast)
spatial criteria (which is really slow)

I would like to run the spatial query only on the filtered layer (which is only one polygon). Is this possible?

SELECT 31H.* FROM bois_corridor_forestier, 31H
WHERE bois_corridor_forestier.ID="1" /*run this first*/
AND st_intersects(bois_corridor_forestier.geometry, 31H.geometry)


texts: Export large area data in Google Earth Engine

texts:
I am trying to average an image collection from 2019-02-01" to "2019-02-28 on Earth Engine platform. My goal is to export the data in Google Drive that I can use it in QGIS for other processing. But several time I got an error message. It says ''Error: Export too large: specified 338272704 pixels (max: 100000000). Specify higher maxPixels value if you intend to export a large area.'
var aoi =
    /* color: #d63000 */
    /* shown: false */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[85.93714952926396, 27.376833254008353],
          [85.93714952926396, 21.187170010173027],
          [89.90596056442021, 21.187170010173027],
          [89.90596056442021, 27.376833254008353]]], null, false);

// get the data
var y2019 = ee.ImageCollection("COPERNICUS/S5P/NRTI/L3_NO2")
              .filterDate("2019-02-01","2019-02-28");

// add layers to map
var nox = y2019.median()
               .select("NO2_column_number_density")
               .clip(aoi)
Map.addLayer(nox,{min:0.00002,max:0.0005,palette:"lightblue,orange,yellow,red,purple"},"NO2_19");
Map.centerObject(aoi, 6)

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: nox,
  description: 'NO2',
  scale: 30,
  region: aoi,
})

You can find the link here


texts: Projecting point in Leaflet from EPSG 900913 (same as 3857) to EPSG 2232 with proj4js returning unexpected results

texts:
I have a leaflet map from which I want to be able to query an API that accepts x y values in EPSG 2232 as inputs. Leaflet doesn't handle point transformations like OpenLayers does so I am trying to use proj4js. I created an epsg2232.js file for the 2232 transform and am using the existing epsg900913.js file that shipped with proj4 (see code below for actual EPSG strings). The general work flow is:

click event returns lat/long in Leaflet default EPSG 3857 (same as 900913)
convert lat/long to EPSG 2232 x/y (see code below).

The code below currently does perform a transform but the numbers are off. For example:

lat = 39.74043 long = -465.01719

transforms to

X = 40674737.46566756 Y = 10890477.404555319.

I expect the transformed x/y to be more like this:

x=3135764.7952012 y=1694936.04430735.

I am unsure what the issue is as this seems like a simple task. It is very possible that I am missing something obvious. Does anyone have any idea what is causing the transform to behave this way?

// formatting lat long click event results to pass to proj4
// not the best method but bare with me
var lat_lng = e.latlng;
var lat_lng = lat_lng.toString();
var lat_lng = lat_lng.substr(7,20);
var lat_lng = new String(lat_lng);
var epsg900913 = new Proj4js.Proj('EPSG:900913');
var epsg2232 = new Proj4js.Proj('EPSG:2232');
var xy = new Proj4js.Point(lat_lng);
Proj4js.transform(epsg900913,epsg2232,xy);

// here is the EPSG string that epsg900913 references.
//Proj4js.defs["EPSG:900913"]= "+title=GoogleMercator +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs";

// here is the EPSG string that epsg2232 references.
// Proj4js.defs["EPSG:2232"] = "+proj=lcc +lat_1=39.75 +lat_2=38.45 +lat_0=37.83333333333334 +lon_0=-105.5 +x_0=914401.8288036576 +y_0=304800.6096012192 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs";


texts:
I need to ask you guys after reading this article.

Any speculations on why some maps showed it?

What can be the cause of inaccuracies in (web)Maps?


texts: Where can I find GIS files that shows the 11 nations of the USA?

texts:
This book describes the cultural boundaries that define America.  At least one person has created what seems to be a GIS style map from this.

Does anyone know where I can locate these files?


texts:
I would like to use QGIS to display solar data from NREL. The file is supposed to provide solar radiation data of Colorado.

However, I have some trouble to select this data in QGIS:
In QGIS, I select "Add Vector Layer" and select the .shp file. QGIS displays the raster for the information (10 km by 10 km data). All squares have the same color. But I can't see the actual data.

I was always under the impression that the shapefile would also contain the geographical data. Is that correct? Do I need to load any other projection file in QGIS to display the solar radiation per latitude and longitude?

How can I display the NREL data correctly with QGIS?


texts: How does a PostGIS query execute from database to the client side?

texts: I am trying to understand how a query is executed from the postgis database to the client. Can someone explain this to me in layman's terms?


texts: Code from 3067 to 3857, not drawing polygon

texts:
Below is a copy of my function, its working without errors. I consloe.log the coords and they are correct, but obviously in 3067 format.

If I remove the code for projection, the 3857 coordinates already saved in database show properly in the polygon, so I am sure that part of the code works. But, I need to project some 3067 coords here also.

I tried a few kinds of projection code from this site but none worked as expected.

I am quite new at mapping code.

  var raster = new TileLayer({
                source: new OSM()
             });

            var coords=this.$store.state.MapsSingle.item.coordinates
            var geojsonObject = {

            'type': 'Feature',
            'geometry': {
            'type': 'Polygon',
            'coordinates': [JSON.parse(coords)]
          }
                      };
        var source = new VectorSource({features: (new GeoJSON({
            dataProjection: 'EPSG:3067',
            // featureProjection: 'EPSG:3857'

//I tried these both ways. exchanging the values for dataProjection and featureProjection

        })).readFeatures(geojsonObject)});

      var vector = new VectorLayer({
        source: source,
        style: new Style({
          fill: new Fill({
            color: 'rgba(255, 255, 255, 0.2)'
          }),
          stroke: new Stroke({
            color: '#ffcc33',
            width: 2
          }),
          image: new CircleStyle({
            radius: 7,
            fill: new Fill({
              color: '#ffcc33'
            })
          })
        })
      });

      var map = new Map({
        layers: [raster, vector],
        target: 'olmap',
        view: new View({
          center: [339042,6673664],
          zoom: 4,
          projection: 'EPSG:3857'
        })
      });
      console.log(map);


texts: Extend “reduceRegions” function application for more the 5000 elements in GEE

texts:
Does anyone knows if is it possible to extend “reduceRegions” funtion application of GEE for more than 5000 elements?

I think it is by default limited to 5000 points/polygons in a feature collections, when I call it as following:

var Xfeatures_sum = ImageX.reduceRegions({
  collection: Xfeatures,
  reducer: ee.Reducer.mean(),
  scale: 30
});


texts:
I'm trying to import data from MS Access but unsuccessful. The commands tab (Customize> customize mode> commands) doesn't show "ArcCatalog", so I'm stuck here being failed to add "Add OLE DB Connection" which is required to import MS Access data file into ArcGIS. I'm using ArcGIS 10.3.1.

Can anyone tell what to do now to import MS Access data into ArcGIS?


texts: Missing grib metadata when using CDO tool to invert Lat axis from South-North to North-South?

texts:
I have a grib file (it is quite big) and I tried to use cdo to copy it to another grib file.
cdo -R copy ws10m_mean24h.grib test.grib

However, in this diff: https://www.diffchecker.com/w25sayVz it shows that the output of CDO on the right hand side misses the important metadata (Ensemble) and dataDate and startStep are altered.
Could someone explain to me what happened and how can I preserve the change from the source grib file: ws10m_mean24h.grib?
What I wanted to do is to invert the Lat axis (from South-North order to North-South order) of ws10m_mean24h.grib to test.grib but all the medata must be kept from ws10m_mean24h.grib.
cdo invertlat ws10m_mean24h.grib test.grib

Thanks,


texts:
I've got a lot of single-band Geotiff files with the GDT_Float32 format and I need to present these over the web, hopefully using a webmap leaflet/mapbox/googlemaps widget where clients can hover over the map and explore the values.

I've been looking at gdal2tiles.py but it resamples everything down to 8-bit png. Converting to 8-bit isn't going to work for me since I need to retain the precision.

Is there a solution out there for creating 32-bit png tiles for a web viewer? Should I be thinking about this another way?


texts: Calculating Cumulative Distance (incorporating slope) of Line Feature Using ArcGIS for Desktop?

texts:
In ArcMap, when creating a line feature class within a geodatabase, the length of each record is automatically derived upon creation. This length is based on a distance measurement for a planar surface. Let us say, however, that the line feature is meant to represent an entity in an area with high topographic relief.

Wouldn't the actual length of that line be different than the result from this planar measurement?
Wouldn't the cumulative length of a line be longer than the planar distance measurement anytime the line does not sit on a planar surface?
How much of a difference is there between the planar and cumulative length measurements?


texts: OpenLayers.Renderer.Symbols custom named by function

texts:
Is there any approach to define custom OpenLayers.Renderer.symbol by a function?

For example, instead of making a symbol for each 3-pointed, 4-pointed, 5-pointed etc. star

OpenLayers.Renderer.symbol.star3 = [ /* coordinates array of 3-pointed star */ ];
OpenLayers.Renderer.symbol.star4 = [ /* coordinates array of 4-pointed star */ ];
OpenLayers.Renderer.symbol.star5 = [ /* coordinates array of 5-pointed star */ ];
...

I want to shorten the code and create a group common symbol defined by a function, like:

OpenLayers.Renderer.symbol.pointedstar(n) = function(n) { /* return array of n-pointed star */ }

and then have it in graphicName style definition:

starstyle = new OpenLayers.Style({ graphicName: "pointedstar(3)", ... });

and even make it depending on feature properties:

starStyle = new OpenLayers.Style(
    { graphicName: "${numofpoints}", ...},
    { context:
         { "numofpoints": /* features contain a property "numberofpoints" to use */
             function(e) { return pointedstar(e.feature.attribute.numberofpoints); }
         }
    });


texts: convert vector (point shape) to raster (tif) using python gdal lib in qgis

texts:
I'm trying to convert a vectorlayer (point- shapefile) to a raster-dataset (tiff) using python gdal-library in qgis.

On Graphical UserInterface I simple choose the menue "raster->conversion->rasterize".
To get what I want, I fill in the file names in tool dialog box and choose "resolution in map units per pixel"
horizonthal: 0.033 and vertical: 0.0164.
That works fine, but when I try to script I orientated on documentation and adapted the code from source:
http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#convert-an-ogr-file-to-a-raster
like this:

# 1. Define pixel_size and NoData value of new raster
NoData_value = -9999
x_res = 0.03333378
y_res = 0.01666641
pixel_size = 1

# 2. Filenames for in- and output
_in = r"C:/Users/.../hoppla.shp"
_out = r"C:/Users/.../hoppla.tif"

# 3. Open Shapefile
source_ds = ogr.Open(_in)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# 4. Create Target - TIFF
_raster = gdal.GetDriverByName('GTiff').Create(_out, x_res, y_res, 1, gdal.GDT_Byte)
_raster.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
_band = _raster.GetRasterBand(1)
_band.SetNoDataValue(NoData_value)

# 5. Rasterize
gdal.RasterizeLayer(_raster, [1], source_layer, burn_values=[0])

It results in error message:

Traceback (most recent call last):
  File "", line 1, in
  File "C:/Users/.../test.py", line 73, in
    _raster = gdal.GetDriverByName('GTiff').Create(_out, x_res, y_res, 1, gdal.GDT_Byte)
  File "C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\osgeo\gdal.py", line 388, in Create
    return _gdal.Driver_Create(self, *args, **kwargs)
TypeError: in method 'Driver_Create', argument 3 of type 'int'

It doesn't accept float values for x_res and y_res, but how can I enter the proper geometry for my dataset?

Additional Info:

QGIS 2.6.1 Brighton on Win7
crs: WGS84 - epsg4326
input shape horizonthal point distance: 0.03333378 units
input shape vertical point distance: 0.01666641 units


texts: Isolating islands within a network (using QNEAT3)?

texts:
I’ve run into a bit of a pickle. I’m attempting to isolate “islands” within a large OSM-derived network. I’m interested in highlighting areas where the network cant connect too. This of course includes actual, physical islands, but it also includes the interest of my research, which is sections of the network that haven’t been properly created and are not connected to the rest of the network.
I’ve been attempting to do this using QGIS’s QNEAT3 plugin. I’ve been creating interpolations of network costs from a point, and attempting to use that to highlight sections that are not being included in the network. This isn't working however, as network cost calculations include areas from the network, and don’t really highlight segments that are completely disconnected from the network. See the following screenshot of a section of a network connected, and the same section cut out to be a “network island.” The contour lines represent the cost of the network, and as you can see, the island section is not included in my network - which is exactly what I want!

What would the easiest way of highlighting sections of a network that do not connect to the rest of the network be? Is there a way of doing it using QNEAT3?


texts: Filtering attributes in a shapefile using a common ID from an Excel file

texts:
I am using QGIS 2.18.11 and I have a shape layer with lines and polylines I need to review; but not all of them, only certain ones. The ones I need to review are marked in a specific field (rev = 0 or rev = 1) in an Excel file (.xlxs). Both data files share the same ID, so I could connect them.

How I would normally do this is performing a join in my shapefile with the Excel file using the ID; then I would add a new field in the shapefile and update it with the "rev" field data; I would then close the excel file so I can keep on working on it in parallel; finally I would filter my shapefile with "rev"=1, and there I have the lines and polylines I need to review.

Is there a faster way? Also a way that does not make me edit the shapefile (by creating that "rev" field)?


texts: Exporting Multiple Data Frames with ArcPy and Data Driven Pages?

texts:
I have a problem in Data Driven Pages
In my Mxd proyect I have 3 Data frames, each frame has 3 layers with their respective legends. I already activated the process of Data Driven Pages to export a single data frame with Python code, however, I need to export the other Data frames with their layers in png or pdf at the same time.

Can I export multiple data frames (with python code) and their respective layers in jpg or pdf?

This the code to export only one data frame:

mxd = arcpy.mapping.MapDocument("CURRENT")
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
  mxd.dataDrivenPages.currentPageID = pageNum
  arcpy.mapping.ExportToPNG(mxd, r:c\\") ##<--- my path

but, I don´t know how export the other data frames using this code.


texts: Leaflet Draw on Mobile - creating a marker

texts:
So I made a simple leaflet map that has the Leaflet Draw plugin.  I am only allowing users to add a marker. Leaflet is supposed to be mobile friendly, but when I open my website in a smartphone (Galaxy s5 chrome browser app), I cannot touch to create a marker on draw:created.  Is there a specific procedure necessary to adapt the draw plugin for mobile touch events?

I looked at https://github.com/tablackmore/Leaflet.draw/tree/touch and tried to redirect my page to a new index which uses the "touch" draw plugin instead of the original plugin but still nothing is "drawing".  This could be because it is not directing.  I am using this script reference http://www.inmotionhosting.com/support/website/redirects/mobile-redirect.  I want to redirect my site to mobile.html rather than index.html.  Regardless of whether my redirect works or not, shouldn't the basic draw plugin work on mobile?

Does anyone have a working example of a leaflet map on a mobile device with drawing capabilities?

For what it's worth here is my project- but it is a novice example and quite messy: http://terpconnect.umd.edu/~melrose/DVA/index.html

PS. This map shows test data only- no real data is displayed.


texts:
I've got a raster with the cut fill tool, its attribute table shows volume information which I'd like to pass to a polygon overlapping the raster.

The polygon  was originally used to clip a raster file so as to extract all depressions of a terrain, later the cut fill tool was used to obtain material volume to fill these depressions. Every thing shows wonderful on the raster cut&fill output but I haven't been able to pass the volume field back to the polygon file as a new field

I've tried to use raster to polygon, extract values to point, some statistical tools, etc, etc with no success at all


texts: Choropleth with classification from multiple columns in QGIS?

texts:
I want to create a small multiple graph, which is something like this or this.

I know that I can create a choropleth by using the styles panel:

Basically, I pick a column attribute, choose the classification mode, and classify. This results in the legend as pictured above. However, this is bin distribution off of one attribute POP1990, in this case the population in each state in 1990. I have attributes POP1991,POP1992,POP1993, and so on. I can individually classify them, and get 4 different bin distributions. But I want to get the bin distribution off of all POPYEAR columns so that they share one legend (like the small multiple example linked above). Is there a way to do that?

Basically, I want to create 4 individual choropleth maps per year, but I want the legend that they all use to be the same. Is there a way to do this? I know that there is this expression dialog:

But I'm not sure how I would use this to accomplish what I want. I saw this about creating small multiples with virtual layers, but I'm still not sure about the legend.

I'm new to QGIS.


texts:
I have bathymetry data for a river in numpy array format. Unfortunately, I do not have the coordinates for each value of this array, but only the central line where the boat passed (x2). I decided to draw a maximum line 100 meters after the boat's path (x1).

To georeference those I used xmin, ymin, xmax, ymax. Here is a part of the Python script
xmin,ymin,xmax,ymax = [track.coords[-1][0],track.coords[-1][1],up.coords[0][0],up.coords[0][1]]
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)

output_raster = gdal.GetDriverByName('GTiff').Create('/content/myraster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
srs = osr.SpatialReference()                 # Establish its coordinate encoding
srs.ImportFromEPSG(4326)                     # This one specifies WGS84 lat long.
output_raster.SetProjection( srs.ExportToWkt() )   # Exports the coordinate system
output_raster.GetRasterBand(1).WriteArray(np_chan1)   # Writes my array to the raster
output_raster.FlushCache()

With that I managed to geo-reference, but the river curve was ignored. Would like to know how to do to consider this curve of the boat.


texts: How to get mean of pixel values within a buffer zone around current mouse position in QGIS

texts:
for a qgis 1.8 python plugin I want to query the mean of all pixel values of a raster layer within a defined buffer around the current mouse position, not for an existing point layer.

For the buffer sth like this should work:

feat.geometry().asPoint()
(3.5176e+06,5.4062e+06)

buff_geom = feat.geometry().buffer(1000, 5)

buff_geom.asPolygon()
[[(3.5186e+06,5.4062e+06), (3.51855e+06,5.4059e+06), (3.51841e+06,5.40562e+06),     (3.51819e+06,5.4054e+06), (3.51791e+06,5.40525e+06), (3.5176e+06,5.4052e+06), (3.51729e+06,5.40525e+06), (3.51701e+06,5.4054e+06),...]]

but now I try to find a solution to query all raster pixels within the buffer zone and calculate the mean, but I don't know how and if this is possible in QGIS 1.8 at all.
All I found is QgsZonalStatistics, but that doesn't work for my task.

Does anybody have a solution for this or a better idea, how to solve this?

Thanks.


texts: Geometry has incorrect latitude-longitudes

texts:
SELECT ST_SRID(geom) AS s FROM zones GROUP BY s;
returns a single record: 4326
SELECT Find_SRID('public', 'zones', 'geom');
returns 4326
SELECT AsText(geom) FROM zones
returns geometry like:

MULTIPOLYGON(((495468.88595 5464400.7263,495468.89 5464401.38005,...)))

Shouldn't my geometry been standard lat-lngs? Something like:

MULTIPOLYGON(((-122.968967906831 49.31253667013,-122.970356138605 49.3125609975662...)))

I tried executing both of these queries, but neither worked:

SELECT UpdateGeometrySRID('zones','geom',4326);
UPDATE zones SET geom  = ST_Transform(geom);

How can I transform the geometry into my desired format?


texts: How to calculate area of a double-line road network in QGIS?

texts:
I have a line vector file that displays a road network. Most roads consist of a double-line indicating the sides of the road (hence the width of individual segments), but some roads are indicated by a single line. Furthermore, many segments are open shapes, so a line to polygon conversion does not "fill in" the paved roadway correctly, but ends up creating random polygons across multiple road segments.
I'd like to calculate the total area of the paved surface of the roads but am unsure how to do this. If I were to calculate the length of the lines and then use an average width to approximate area, the length would be incorrect because most roads are represented by two lines - but some roads only have one line - so I can't simply divide the total length by two. There is also a lot of variability in the width of the roads.
Is there a way either to convert the lines to polygons - but respecting the "pavement" of the roads? Or another method altogether? I attached a screenshot to illustrate (the buildings are a separate polygon layer).


texts: Reading 3D polygon WKT from MapGuide with Openlayers


texts:
I have modified 20+ shortcut keyins using the QGIS 2.18 Settings > Configure Shortcuts dialog. I save them to an XML file using the Save... button, and the XML file that is generated looks fine (that is, all the modified keyins are displayed).

However, when I load and configure another PC with QGIS 2.18.4, the shortcuts XML I have generated doesn't load on this PC. I click on the Load... button and specify the modified shortcuts XML file, but no changes at all appear in the shortcuts list.

For that matter, I can remove shortcuts for several functions on the originating 2.18 PC, but attempting to reload the saved shortcuts XML file has no effect on this system either. Is anyone else experiencing this problem?


texts:
I am digging into clustering strategies with Open Layers 2. Looking at the example OpenLayers2 Cluster Strategy Threshold, I was able to get it working. Unfortunately I am working with OSM, and the example is not working as expected by replacing the appropriate line in the code example:

 /*var base = new OpenLayers.Layer.WMS("OpenLayers WMS",
                 "http://vmap0.tiles.osgeo.org/wms/vmap0",
                {layers: 'basic'}
            );
 */
 var base = new OpenLayers.Layers.OSM();

It results in looking at only one circle instead of many. Only when I start zooming in, more and more circles do appear. Why is it different when I am using OSM? Which settings have to be taken care of?


texts: Transform wrongly transformed data with .gsb file

texts: I have data which originally came in GK coordinates (EPSG:31467) as an DXF. This data was transformed using the QGIS standard .gsb file for Germany and then got digitalized (basically draw the lines etc into a new layer and fill in attributes which are only visible in the DXF). Now I just came to realize that the transformation was done using the Germany-wide grid, but it actualy should have been transformed using the grid of the state Baden-Würtemberg (BWTA2017.gsb). How can I transform my data which I have already digitized based on the wrongly transformed DXF? Is there any posibility? There is a difference of ~5cm when I load the original DXF into QGIS and reproject it on the fly with the BWTA2017.gsb compared to the German-wide grid.


texts: Calculating AICc value from Maxent output

texts:
I want to calculate the AICc value for my Maxent models. On a little search, I have found that calc.aicc function if implemented in R, can help in finding out the AICc value. But, I am not sure, what should be the input parameters. The link of the function is: https://rdrr.io/cran/ENMeval/man/calc.aicc.html
calc.aicc(nparam, occ, predictive.maps)

get.params(model)

I know that occ is the occurence records and predictive maps might be referring to the ASCII file generated by the Maxent but I am not sure about the other parameters, like what is meant by model here etc.


texts: How to iterate through all fields, and all records in each field with GDAL/OGR

texts:
I'm trying to iterate over each field in a shapefile. Then iterate over each record in that field.

First I build a list of all the fields, then iterate through that list. While I'm going through the list, I loop through all the features via feat = feature.GetField(x) where x is the current field.

This works for the first field. It prints the first field, then all of its records. Then it proceeds to just loop through the rest of the fields, not the records. It's like it just runs the records for loop once. I have no idea why. I think it should work.

from osgeo import ogr

shapefile = r"E:\metadata_2\data\test.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
layerDefinition = layer.GetLayerDefn()

layer_list = []

def getfields():

    for i in range(layerDefinition.GetFieldCount()):
        layer_name = layerDefinition.GetFieldDefn(i).GetName()
        layer_list.append(layer_name)

    print layer_list

    for x in layer_list:
        print 'This is the layer {}'.format(x)
        for feature in layer:
        feat = feature.GetField(x)
        print feat

getfields()


texts:
I want to extract trees from a given normalized Digital Surface Model (nDSM) in SAGA GIS. My workflow is the following:

Applying a Gaussian Filter to the nDSM with search radius with standard deviation=1, search mode=circle and search radius=5m
Reclassifying grid values to exclude vegetation below 5m
Applying an inverse watershed segementation algorithm (searching maxima; not joining the elements)

Somehow, my result yields too many trees a some site, too few at others. Some "real" trees even show up two seed points.

Here's a comparison my found trees to the orthophoto (dark green=trees or pastures. light green=tree shadows/bare):

Besides a little shift due to georeference inaccuracies: Any ideas how to improve the workflow? I already tried to vary the parameters (resampling the 1m resoluton to 10cm, gaussian radius 2-5m, reclassifying threshold=15m (as to extract tree crowns only), joining the segments based on seed to seed difference...)


texts: Placing labels with leader lines at periphery of map?

texts: I have a set of points I want to label with leader lines. I know how to set a defined offset for the leader lines but how do I get Maplex to place the text outside of the map and adjust the leader lines as necessary?


texts: Possible to import current version (v2) of the SQL Server Spatial binary format with QGIS/GDAL?

texts:
I'm running QGIS 3.4.3 on a Windows 10 machine and when loading vector layers from our MS SQL Server version 12 everything seems to be working great. No problems connecting or visualizing the data. But when looking closer at a layer with property lines I realized some are missing. When inspecting the same mssql layer/table with FME Data Inspector the missing lines are visible.

Googleing the issue eventually led me to this bug report/feature request. Summarized it says that QGIS (or acctually GDAL which QGIS relies on for mssql-support) simply cant show SQL Server Spatial binary format version 2 meaning that six year (2012) old changes in my database just can't be visualized!? Version 2 geometries are simply not recognized?

This issue and this pull from Github indicates that this problem was solved with the 3.8.0 update?

But when updating to 3.8.0 I get this error instead and now I can't connect to the server:

[Microsoft][ODBC SQL Server Driver][SQL Server]Login failed. The login is from an untrusted domain and connot be used with Windows authentication.[Microsoft][ODBC Driver Manager]Drivrutinens SQLSetConnectAttr misslyckades. QODBC#: Unable to connect

Am I just being stupid or is there really a way to read SQL Server Spatial binary format v2 with QGIS/Gdal? Has anybody had any success with this? Work-arounds are more then welcome.

------ Additional information by request -------

SQL Server Spatial binary format v2:
I can't access a file to send over, could below documentation do?https://docs.microsoft.com/en-us/openspecs/sql_server_protocols/ms-ssclrt/77460aa9-8c2f-4449-a65e-1d649ebd77fa
https://sqlprotocoldoc.blob.core.windows.net/productionsqlarchives/MS-SSCLRT/%5bMS-SSCLRT%5d.pdf

https://docs.microsoft.com/en-us/openspecs/sql_server_protocols/ms-ssclrt/77460aa9-8c2f-4449-a65e-1d649ebd77fa

MS SQL Server:
Product: Microsoft SQL Server Standard (64 bit)
Operating system: Microsoft Windows NT x64
Platform: NT x64
Version: 12.0.5571.0
Memory: 55282 MB
Processors: 16


texts:
I am very new to Quantum GIS, I want to create a digital elevation map a particular region, I would very much appreciate if some one gives me a step by step procedure.

Thank you

edit by Martin: I have the same question, but I want to create DEM from scan map (the territory of interest is aprox. 28 square km). Work flow: First open the map in an image editor (GIMP), select by colors with some tolerance and select the black color. Invert selection and delete non-black stuff. Save as Tiff WITH NO COMPRESSION.
And then what is next step in QGIS with Grass Plugin? How can I get the contour values to the Tiff (my big problem)? After all is no problem make DEM from contours, I guess.
Sorry for this "stupid" questions and thanks for all your replyings.
Regards Martin


texts: Aggregating Polygons To Meet Privacy Requirements

texts:
I have a point feature class representing office locations of all employers in a certain industry. The feature class has an attribute for storing the number of employees working at each office. Someone has requested to use this data, spatially-joined to the smallest geographic unit possible — Census Blocks, in this case. However, a privacy agreement prevents the data's release as-is. Instead, it must be suppressed to meet two criteria:

Any polygon must contain at least 3 employers (points);
No more than 80% of the total employment within a polygon may be by a single employer.

I have successfully written a script that spatially joins the points to Census Blocks, keeping the sum and maximum employment in each. Each one that does not meet the suppression criteria is flagged. (Polygons containing no points are not flagged, since there is no data to suppress.) I then check each Block Group to see whether any flagged Blocks are contained within it. Block Groups containing only un-flagged Blocks are then replaced with the Blocks. The resulting features class is then checked against the suppression criteria, to check whether the Block Groups have adequately suppressed the data.

The same process is repeated for Tracts, leaving me with a dataset consisting of Tracts (some flagged and some not), Block Groups and Blocks (all un-flagged). The next progression in the geographic hierarchy, however, is the county, which is of no use to the person requesting this data.

My question, then, is this: Are there any commonly accepted methods of aggregating polygons into as many groups as possible, so that all meet some minimum criteria?

Here are some rules that I would like to apply to the aggregation:

Whenever possible, flagged Tracts should only be aggregated with other flagged Tracts;
For the flagged Tracts that are not contiguous with any others (or isolated groupings that still do not meet the criteria), they can be joined with Tracts that already meet to criteria, although there may be Tracts with no employers in between them that will also need to be included.
I would like to keep county boundaries intact unless absolutely impossible (and I anticipate doing this by separating the inputs features into their respective counties before processing them).
The solution must be in Python, with the use of ArcGIS tools or open-source Python libraries.

Ideally, someone can point me to an existing means of implementing this aggregation. If not, I am happy to code the algorithm myself, although a list of specific steps/tools would be much appreciated. The problem strikes me as a special case of redistricting (with dis-contiguous polygons), and to this end I have looked into using PySAL's regionalization algorithms, although it is not clear to me how to check the maximum employer's percentage of total employees using these.


texts: QGis 3.4.4 Create Atlas - one geometry per page

texts: I have a shapefile with 50 Polygons. I want to create one PDF file in which I see only one polygon per page. I think I have to use the Atlas function. But I don't know how to Setup the filter.


texts: ArcGIS messing with feature geometry in Editing mode

texts:
I am currently working on project, where I need to digitise a geomorphological map. I am working in ArcMap 10.5.1. I have an issue with creating new features. When I am digitising an object from the map sometimes, a strange error occurs with the geometry. It looks like after drawing some vertices, something's going wrong, and it seems like ArcMap on its own finished part of the sketch and closed the part of the polygon in very strange way. Here's the example:

When I was trying to fix the polygon manually with Edit vertices tool, it didn't worked well, because those long random lines are intersecting the part I have already created and geometry gets really complicated, so it's usually impossible to fix it manually.
What's more interesting, during editing the veritices I have realised that the beginning vertex (which is indicated in ArcMap as red) have changed it's location, and now different vertex is indicated that way. Below here is example - I have indicated with yellow arrow the original position of begining vertex.

Sometimes, ArcMap is not "finishing part on its own", but those strange lines appear. Also, what's important it doesn't appear with every new feature I am creating, but it is often repeating in the same places of the map (e.g. I was trying to draw a new feature in the same place, when last time this error occured, and it occured again).
I have been looking for similar problem, but I haven't found any. I also have tried the classical solution with reinstalling the software but it didn't work. Any ideas what can cause such problem?


texts: Distorted h3 hexagons using h3js with Mapbox

texts:
When a user clicks on a Mapbox streets-v11 map, I'm attempting to:

Take the clicked point, and identify the h3 hexagon (at resolution 10) around the point
Convert the hexagon coordinates into the polygon hexPolygon
Highlight building features (from streets-v11/composite source) within hexPolygon

Currently, this functionality is working as expected, except the h3 hexagon coordinates are forming a narrow distorted hexagon.  Below is an image of the distorted hexagons being created.
Here is a jsbin with my code.  To reproduce, can zoom in on the map until you see buildings render, and click on the map.  At that point, a narrow hexagon of buildings should appear highlighted.
(You can then copy the hexPolygon object from the console and paste it into a tool like geojson.io to view the polygon that was created)
And here is the function:
map.on("click", function (e) {

  // Convert a lat/lng point to a hexagon index at resolution 10
  const h3Index = h3.geoToH3(e.lngLat.lng, e.lngLat.lat, 10);

  // Get the center of the hexagon
  const hexCenterCoordinates = h3.h3ToGeo(h3Index);

  // Get the vertices of the hexagon
  const hexBoundary = h3.h3ToGeoBoundary(h3Index);

  hexBoundary.push(hexBoundary[0]);

  console.log(h3Index);
  console.log(hexCenterCoordinates);
  console.log(hexBoundary);

  const hexPolygon = {
    type: 'Feature',
    properties: {
      paint: {
        "line-color": "red"
      }
    },
    geometry: {
      type: 'Polygon',
      coordinates: [hexBoundary]
    }
  };

  var features = map.querySourceFeatures("composite", {
    sourceLayer: "building"
  });

  console.log("features:");
  console.log(features);

  const collector = [];

  console.log('hexPolygon:');
  console.log(hexPolygon);

  for (let i = 0; i < features.length; i++) {

    const feature = features[i]

    if (feature.geometry.type === "Polygon" && turf.booleanWithin(feature, hexPolygon)) {
      collector.push(features[i]);
    }
  }

  console.log("collector:");
  console.log(collector);

  const buildingsHighlightedLayer = map.getSource("buildings");
  buildingsHighlightedLayer.setData({
    type: "FeatureCollection",
    features: collector
  });
});


texts:
I have a osm file,i convert it to pbf file by osmosis.but it not convert to .map file.i use this in cmd:

osmosis --rb
file=city.osm.pbf -mw file=map.map bbox=34.3635451292,30.0247495993,57.513131979
3,53.2954765066

get it ERROR:

SEVERE: Execution aborted.
org.openstreetmap.osmosis.core.OsmosisRuntimeException: Argument bbox for task 1
-rb was not recognised.

what is this?!

My Osm File Bound Tag:

<bounds minlat="30.0247495993" minlon="53.2954765066" maxlat="34.3635451292" maxlon="57.5131319793"/>


texts: GeoPandas GeoDataFrame plot statistics - how?

texts:
The plot method of a GeoPandas GeoDataFrame overrides the DataFrame parent method to plot the geometries of the associated GeoSeries. This is very convenient, but I'd also like to use the easy ways Pandas provides to plot other properties of my dataset.

However, when I now try to make a scatterplot, or any other kind of plot, I get an error message:

mygeodataframe.plot(kind='scatter', x='somecolumn', 'y'='someothercolumn'...)
*** TypeError: plot_dataframe() got an unexpected keyword argument 'kind'

How do I get my plotting methods back? Converting a copy of the GeoDataFrame (minus the 'geometry' Series) into regular DataFrame helps, but it's clumsy.

I'm using GeoPandas version 0.1.1 and Pandas version 0.17.1.


texts:
I would like to build a web-based application using leaflet which will have a flat plane and a simple X,Y coordinate system. Is there an SRID I can use for this? I will be generating my own basemaps based on the coordinates.

If it helps you conceptualize what I am doing just consider the game Zork where you move in an X,Y grid. It is similar except I would like to use Leaflet to create this to allow the user to pan, zoom and load basemaps for the different Z-Indexes.


texts: How to interpret the legend using the Flannery method with the size assistant in QGIS?

texts:
I've tried to use the data defined override with the Size Assistant as follows
Scale Method Flannery.
I'd like to know how this method is working and how can we read the legend generated automatically by QGIS (see figure ):

Data defined override
Active: yes   (ctrl|right-click toggles)
Expected input:
double [>= 0.0]
Valid input types:
string, int, double
Current definition (expression):
coalesce(scale_exp("fieldA", 2.11, 3.52, 1, 15, 1), 0)


texts:
It seems that when geometries are removed from an IGraphicTracker instance, the geometries remain in memory until IAoInitialize.Shutdown() is called. I have tried removing the geometries by calling the IGraphicTracker.RemoveAll() method and by storing the IDs returned by IGraphicTracker.Add() method and calling the IGraphicTracker.Remove() method for each ID.

Is there something that I'm missing in the example provided below?

I'm currently using ArcGIS 10.1 and Visual Studio 10 targeting .NET Framework 3.5.

An example can be found here: http://forums.arcgis.com/attachment.php?attachmentid=20140&d=1356115382

This problem was originally noted on the ESRI Dev Forum: http://forums.arcgis.com/threads/72516-Memory-Leak-with-Graphic-Tracker


texts: Double output with getfeaureinfo request with Mapserver and Openlayers

texts:
I am trying to make an popup InfoTool in OpenLayers with a Mapserver file. Everything is working fine now. But when I click on the map I gives a double output. So I get for instance: cityname, citycode, cityname, citycode

I cannot find the problem. I think it has something to do with the mapfile.

MAP # start ************************************

NAME "bag"
EXTENT 6180 300854 287160 629230
STATUS ON
DEBUG ON
UNITS METERS
CONFIG "MS_ERRORFILE" "mapserv.log"
CONFIG "PROJ_LIB" "c:\proj\nad\"
TRANSPARENT ON
SYMBOLSET "symbols\agnsymbols.sym"
FONTSET "fonts\fonts.list"

PROJECTION
    "init=epsg:28992"
    "ellps=bessel"
END

WEB
    IMAGEPATH "\tmp"
    IMAGEURL "\tmp\\"
    #TEMPLATE "template\woonplaats.html"
    METADATA
        WMS_TITLE "Geoserver.nl BAG service"
        WMS_ABSTRACT "Geoserver.nl BAG service"
        WMS_ACCESSCONSTRAINTS "none"
        WMS_ONLINERESOURCE "http://geo.nl/bag/mapserv.cgi?map=test_stefan.map&service=wms"
        WMS_SRS "EPSG:28992"
        #WMS_ENABLE_REQUEST "*"
        #OWS_SCHEMAS_LOCATION "http://ogc.dmsolutions.ca"
        WMS_FEATURE_INFO_MIME_TYPE "text/html"
    END
END

OUTPUTFORMAT
    NAME png
    DRIVER AGG/PNG
    MIMETYPE image/png
    IMAGEMODE RGB
END

# Start of LAYER DEFINITIONS ---
# Woonplaats layer ---
LAYER
    NAME "Woonplaats"
    CONNECTION "user=user password=password dbname=dbname host=localhost port=5432"
    CONNECTIONTYPE postgis
    DATA 'geometrie from agn.woonplaats using unique wnpcode using srid=28992'
    METADATA
        WMS_TITLE "Woonplaats"
        WMS_ABSTRACT "Woonplaats"
        #GML_INCLUDE_ITEMS "all"
        #WMS_INCLUDE_ITEMS "all"
        #OWS_INCLUDE_ITEMS "all"
        RESULT_FIELD "woonplaats"
    END
    PROJECTION
        "init=epsg:28992"
    END
    EXTENT 6180 300854 287160 629230
    STATUS ON
    DUMP TRUE
    TYPE POLYGON
    UNITS METERS
    SIZEUNITS PIXELS
    TOLERANCE 0
    TOLERANCEUNITS PIXELS
    TEMPLATE "template\woonplaats.html"
    CLASS
        MAXSCALEDENOM 200000
        STYLE
            OUTLINECOLOR 61 61 61
            WIDTH 4
        END
        STYLE
            OUTLINECOLOR 161 161 161
            WIDTH 2
        END
    END
    CLASS
        MINSCALEDENOM 200000
        STYLE
            OUTLINECOLOR 61 61 61
            WIDTH 1
        END
    END
END

END # End of Mapfile

The template file is as follows:

       Woonplaats
       [item name=woonplaats format=$value escape=none]
       [item name=wnpcode format=$value escape=none]

Hope someone can give me some more insight


texts: Applying template (layer file) to raster in ArcGIS ModelBuilder?

texts:
I am trying to apply a template layer, .lyr file, to a raster created as an output of raster calculator in ModelBuilder.

However I cannot set this as the input feature to the tool "Apply Symbology from layer".

See below screenshot:

What is the error in my approach?


texts: ST_Difference in PostGIS : problem with geometry result (Polygon)

texts:
I'm experimenting a result issue using the ST_Difference spatial function in PostGIS.
Input :
Feature 1 and 2 overlaps each other. Both contains 5 vertices, like a rectangle (4 plus an additional one to close the polygon).
I want to get the geometry difference between feature 2 and feature 1. (st_difference(geom(feature2),geom(feature1))

Result :
After executing the function, the result of geom(feature2)-geom(feature1), in green, contains unexpected residual vertices (highlighted), inside feature 1 :

Result only :

Anyone knows how to avoid this issue, or how to get a "cleaned" difference of 2 input polygons?

The input data is Polygon, the result (after st_difference) is multipolygon. I  get the same result with Polygon geometry for the result.
Input :

Result :

Edit 2 :
Input polygon geometry :
Feature 1 :
"POLYGON((363011.683640322 5784798.67947383,363061.369728361 5784794.01684973,363060.058365334 5784774.34640432,363009.64374228 5784780.75751245,363011.683640322 5784798.67947383))"
Feature 2 :
"POLYGON((363079.583103742 5784786.58579258,363079.583103742 5784777.55195839,363053.501550197 5784780.32039144,363054.230085212 5784789.06281163,363079.583103742 5784786.58579258))"
Result of st_difference(feature2(geom),feature1(geom)) :
MULTIPOLYGON(((363079.583103742 5784786.58579258,363079.583103742 5784777.55195839,363053.501550197 5784780.32039144,363060.407760381 5784779.58733003,363060.995393925 5784788.40183319,363054.230085212 5784789.06281163,363079.583103742 5784786.58579258)))
As shown on the screenshot, the result contains 7 vertices (instead of 5 as expected).
CRS : 28355


texts:
I am attempting to run the SAVI function (found within the Spatial Analyst module) in ArcPy. However, this function is not showing up for me regardless of there being multiple ArcGIS pages with code samples. (https://pro.arcgis.com/en/pro-app/arcpy/spatial-analyst/savi.htm)
Is this function only available for ArcGIS Pro users?
I have a Spatial Analyst license and I am using ArcGIS Desktop. Is this a licensing issue? Why can't I execute the SAVI function in my code?
SAVI_raster = arcpy.sa.SAVI(Raster,5,4,L_Value)

Error message:

AttributeError: 'module' object has no attribute 'SAVI'


texts: QGIS driving instructions predetermined route

texts:
I have a single line layer (a predetermined route), extracted vertices with distance, angle and added streetnames.
I now want to add instructions (automated), like a regular sat nav would do.
For  starters  with the simple 'turn right', 'turn left' .
My goal is to export the file to PostgreSQL, and use a mobile device to drive the route. (using OpenLayers)
Any ideas and/or experience?


texts: Understanding the performance boost behind type casting in PostGIS

texts:
I am running a PostGIS server (version 2.3). I have a table with the following structure:

- index: BIGINT
- geometry: LINESTRING

I was playing with some of my spatial queries (ST_Intersects in this case) to understand the performance better. I had a query that looked like

SELECT * FROM my_table WHERE ST_Intersects('SRID=4326;POINT(x, y)', my_table.geometry)

As you can tell, table structure and query itself are all simple. Note that geometry column has only LINESTRINGs that follows the same coordinate system (referred with SRID 4326).

As soon as I updated my query with a type cast (SELECT * FROM my_table WHERE ST_Intersects('SRID=4326;POINT(x, y)'::geometry, my_table.geometry)), the performance improved by 10-20x.

Is this happening because PostGIS doesn't need to find out if input geometry is geometry or geography? Not sure, because the second argument of ST_Intersects is guaranteed to be a geometry, because it was created with the following query: SELECT AddGeometryColumn('my_table', 'geometry', 4326, 'LINESTRING', 2);.

Does anybody have any insight on what would be reason for the performance boost with the type cast?


texts: Content for Introduction to GIS course using QGIS?

texts:
There is already a Q&A here at GIS SE on Where to find QGIS tutorials and web resources? but this is not a duplicate of that because it is about Instructor-led rather than On-line Training.

For teaching an Intro to GIS class using QGIS, what could be its content?

As an example, the class could be for an Masters in Public Administration program, where students:

often look for jobs as government analysts and managers, or non-profit executive directors
need to be able to make a simple map that can go in a document
need to be able to talk to technical GIS staff intelligently
will probably only ever take one course in GIS
probably won't be looking for work as GIS techs

Such a course could have these requirements:

20 hours material to be completed over a weekend
includes a couple of hours demonstration on QGIS to ESRI parallels
to be taught using QGIS

My basic text will actually be "Mapping It Out" (Monmonnier, 1996 ??), which they will read before class and write a short page on, since it is most important they learn how to make an argument from a map as part of a bigger decision making process.  I will assign the QGIS tutorial to be skimmed, and have as optional Glaser's book on QGIS and old editions of Dent's Cartography.

I think I will do a demonstration with questions (~2 hours), a little bit of lecture on basic concepts (~3 hours, distributed over the weekend), a lot of practice time (~12 hours), and then student presentations and discussion (~3 hours).

The demo and student maps will start with some addresses, make a locator map of the points, edit them a little for practice, do something with the field calculator, figure out the local demographic characteristics of the points by which tract they are in, calculate some point-in-polygon rates based on local population and point/ event counts, and finish by embedding the PNG in a document and writing a paragraph that refers to the map.


texts: Exporting Ui.chart in Google Earth Engine?

texts:
I have a project where the user press on a district and then the app shows the temperature a s time series chart. I want to optimize all that calculation and save each chart instead of making the whole computation each time.

For that I need to export ui.chart like this :

var tempchart = ui.Chart.image.seriesByRegion(CollectioninCelsius,
                       district,
                       ee.Reducer.mean(),
                       'Celcius',
                        1000,
                       'system:time_start',
                       'PROJECT').setOptions(title);


texts:
It may be a bit off-topic. I am eager to know the industry specific importance of Spatial Databases.

Is Spatial Database a very important part of any GIS system? Or will we use other technologies to store and retrieve spatial data in the near future?

I want to know these things because soon I am going to pursue MS on Spatial Database Optimization.


texts: an arcgisscripting method for sanitizing a filename?

texts: Is there an arcgisscripting method for sanitizing a filename to get rid of characters which Arcgis barfs on, the meek hyphen - for example?


texts: What projection / coordinate system are used for Maps.me?

texts:
I need to import data into Maps.me but since I assume there is not OTF reprojection, I need to know the CRS I must use.

Any ideas or suggestions?


texts: How to visualize measure values along linestrings in QGIS?

texts:
I have a series of LineStringM geometries. The measure value represents slope (or stream gradient) at each point in the geometry. I love QGIS's robust styling options, but I can't find a way to render each measure value along my linestrings. Ideally, I would have blue for slow gradients, and red for very steep gradients.

To help understand what I'm trying to accomplish, let's suppose we've used Linear Referencing to encode the length of the LinestringM thus far into each point's measure value.  At the top, you would have 1.0 and at the bottom, you would have 0.0.  I've added an image what what I'm hoping to accomplish below. I made this in Paint.NET, so it's just a mock-up.

direct link to image

How would you go about accomplishing this in QGIS?


texts: How to use ArcObjects to limit classes visible in legend?

texts:
Is there a method use ArcObjects for .NET to manipulate the "Only show classes that are visible in the current map extent" option for legend items?

The GUI for this functionality is accessible via the legend properties:

I have read through the ESRI documentation for the interfaces supported by Legend, and ILegendItem3 but have not found anything relevant.


texts: Problem in reading mosaic_dataset path from a csv file

texts:
Based on this How to batch Check geometry in ArcGIS for a list of datasets from a CSV? code, i'm trying to extract the mosaic dataset path's within multiple folders. If i directly provide the mosaic_dataset path instead of fcs, the code works fine. But don't know why it's not reading the path from csv file. This is the path in the csv file(row 0) D:\QAQC\Test.gdb\Aus_Geology.gdb\250K_Geology. If i copy and replace fcs with exact same thing, the code works.

import csv
import arcpy
csvpath = r"D:\QAQC\Open2.csv"
with open(csvpath, "r") as csvfile:
    fcs = [r[0] for r in csv.reader(csvfile)]
arcpy.ExportMosaicDatasetPaths_management(fcs,"D:\Test\Paths.dbf")

This is the error i am getting

File "", line 6, in
File "c:\program files\arcgis\desktop10.3\arcpy\arcpy\management.py", line 11099, in ExportMosaicDatasetPaths raise e
RuntimeError: Object: Error in executing tool


texts:
I am having difficulty writing a nested postgis query. I have a single table from which I need to create a view in PostgreSQL.

I have created the view with the following command

create view rrop08.rrop08_lws_f AS
select * from rrop08.rrop08_gk3_polygon
where planznr='04.01' OR planznr ='04.02';

The problem here is that the table contains another field called "bez". I need to carry out the first query but also with bez=0; so that none of the 04.01 or 04.02 rows have a bez value other than 0.

I have tried this...

create view rrop08.rrop08_lws_f AS
select * from rrop08.rrop08_gk3_polygon
where planznr='04.01' OR planznr ='04.02' AND bez=0;

but this retrieves all rows with bez=0 and not just the 04.01 or 04.02 ones. I need to first select 04.01 or 04.02 and then on this query get those with bez=0...ie.nested.

I´m a bit lost on how to do this.

THanks for any help,

Rob


texts: All of my scripts in GEE code editor have disappeared

texts:
I signed into GEE code editor today to find that all of the scripts I've written and saved are now gone. It just says "No accessible repositories" under owner/writer/reader drop down panel. All of my assets are still there though. I do use a government email address to sign in if that matters.

Is there anything I can try to get them back?


texts: Projection problems when importing KML from Google Earth to QGIS

texts:
I want to open several points that I exported from Google Earth as a KML in QGIS for further use. I can import the data and the points are displayed but something must be wrong with the projection as the points are plotted somewhere completely different.

When I open the KML separately the projection of the file is in WGS84 EPSG:4326 (which was also the setting of the original project). I enabled projection on the fly and also tried exporting the KML as a shapefile in different projections but it didn't work.

Is it maybe due to the fact the data are formatted in degrees / minutes and how can I adjust that?

I use QGIS 2.0.1 on a Mac.