texts: Interpolate points against angled polygons

texts:
I have a shapefile, made up of many individual polygons that are essentially arranged in a pixel structure (this file: https://data.gov.uk/dataset/vessel-density-grid-2014). I am interested in the Annual_Ave field data.

I also have a regular array of points that I wish to interpolate against this data. The problem that I have is that the vector data is at a different angle to my array of points (my points are parallel to lines of constant latitude and longitude). The first image roughly illustrates what I mean, the second image is a close up of the actual data to illustrate the problem. I want to be able to perform this interpolation in Python, using GDAL/OGR, and am happy to store the vector data in a different format to facilitate this.

I tried converting the vector polygon layer to a raster, using the ArcGIS polygon to raster function, in the hope that I could extract the latitude and longitude points of each pixel and use this to interpolate my own data. However, when I import the resulting raster into Python using GDAL, the rotational information does not seem to be present so the geotransform is incorrect.

Example (raster file created with ArcGIS Polygon to raster, Cell Size=0.004):

import gdal
in_ds = gdal.Open(<path to raster file>)
in_ds.GetGeoTransform()
Out[39]:
(-10.036726874521381, 0.004, 0.0, 63.932379486900615, 0.0,
-0.003999999999999999)

As you can see, the rotation field is 0 even though the polygons are clearly rotated as can be seen in the file. This means that the latitude and longitude values that can be generated from the geotransform data do not match up with those for the pixels in the file.

The other thing that I considered was loading all of the polygons into Python, getting the latitude and longitudes of all the points from the polygon centres, ordering them and interpolating my data points against this. The issue I think with this is that it will be pretty slow, as it is a large vector dataset containing some 280,000 polygons.

Ideally, I think that it would be preferable to save the vector layer as a rotated, interpolated raster as it is essentially raster data that they have stored as a vector for convenience.


texts: Getting the average of maximum values from raster (tree top segementation)

texts:
I've got raster with elevations of trees in a forest. I would like to find out the average tree height of all trees in a given area.

Because of the geometry of a tree top (which is like a cone and just the top is my real tree height), common statistic tools (Zonal Statistics or v.rast.stats) doesn't work for me. I need the average of the maximum values.

The raster:

All trees and vegetation below 4m have no data. Even all values above 40m. The yellow area in the north-west seems to have mostly 40m heigh trees, but the statistics (median/mean) say something about 28m.

Any ideas, how to approach this task using QGIS, SAGA, GRASS?


texts: VBscript to change label only for certain values

texts:
I want labels displayed for two fields: [ID_Messpkt] and [Pb_corr], both numeric type
I want to change the label for certain values of one of the fields: [Pb_corr]

I have two scipt that work individually but I don't know how to group it in one script. Does anyone know how to do that please?

Function FindLabel ( [ID_Messpkt], [Pb_corr] )
   FindLabel = [ID_Messpkt]& vbNewLine &[Pb_corr]

 Function FindLabel ([Pb_corr] )
   if ([Pb_corr] = 49) then

     FindLabel = "<50"

   end if
 End Function

Found the easy solution, sorry

```
Function FindLabel ( [ID_Messpkt], [Pb_corr] )
   FindLabel = [ID_Messpkt]& vbNewLine &[Pb_corr]

   if ([Pb_corr] = 49) then

     FindLabel = "<50"

   end if
 End Function


texts: In MapBox Studio, how to create a centred label of a polygon created from shapefile?

texts:
To get the undesired result below, I:

Uploaded a zipped shapefile to MapBox Studio.
Using the shape, created a layer of polygons to display area boundaries.
Using the same shape, created a text layer to label areas with their names.

Labels appear on the lines:

How can I make labels appear in the centre of the polygons?


texts:
I am using QGIS and I have an Excel file that has fields with direction and it fails to geocode these addresses.

Examples:

address                apartment     township
1 calle 2 avenida      Guatemala     Escuintla
4 calle 5 avenida      Guatemala     Santa rosa


texts:
I am attempting to extend a series of polyline features to an adjacent polyline feature. I need the lines to intersect, but not extend beyond the 'first hit' of the adjacent polyline feature and maintain their original bearing. As such, distance bearing to line using a set distance and the edgematch tools aren't suitable. The extend lines to layer tool (Extending batch of lines to another object using ArcGIS Desktop?) just crashes/hangs. And ETGeoWizards/Geotools doesn't seem to have the necessary tools. Any other suggestions. I would like to create a point layer at the line intersects.
I'm using ArcMap 10.5.1 with advanced license. (Image: Yellow lines = line segments to extend, pink = adjacent shoreline. Red dashed line = what I would like to happen)


texts: Adding Leaflet legend based on symbology?

texts: I have a question related to legend. I have a series of layers as the picture below. I want to add legend (symbology) base on categories on that map. For example, if I select layer land use, the symbology (color and categories) will appear on that map.


texts:
I have two lists of coordinates in lat/long format which are showing the same points but in different coordinate systems.

I need to find the corresponding points in two datasets using distance but as they are in different coordinate systems, it is impossible.

On the other hand, the coordinate systems are not known so that I can use a geographic transformation.

Can I use Affine or projective transformation by extracting few tie points between two datasets?


texts:
I'm having trouble selecting individual the corresponding polygon in my parcel data based on the location of a property. It seems obvious, but what am I missing?

I have a list of geocoded properties labeled as "active". Then I also have polygon parcel data. I'd like to iterate through the list of active properties and determine which parcel and the accompanying data that corresponds to.

I am trying to use Select by Location, however it is not returning any values. Any thoughts?

Model builder that iterates through the active properties, then Selects the correct polygon based on the Selected feature from active in the iteration

Image of the parcels and the active properties


texts:
I am using the latest version of QGIS (Lyon) on MAC and have set the default Python version 2.7. Moreover, I have installed all the required python modules and GDAL.

It was working before, but now, when I click on the "Browse" button in

Vector-> Geoprocessing->Union

it returns the following error (actually, it returns this error whenever I click on any "browse" button in QGIS):

An error has occured while executing Python code:

  TypeError: QgsEncodingFileDialog(QWidget parent=None, QString caption=QString(), QString directory=QString(), QString filter=QString(), QString encoding=QString()): argument 3 has unexpected type 'QPyNullVariant'

  Traceback (most recent call last):
   File "/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools/doGeoprocessing.py", line 125, in outFile
     (self.shapefileName, self.encoding) = ftools_utils.saveDialog(self)
   File "/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools/ftools_utils.py", line 327, in saveDialog
     fileDialog = QgsEncodingFileDialog(parent, QCoreApplication.translate("fTools", "Save output shapefile"), dirName, filtering, encode)
  TypeError: QgsEncodingFileDialog(QWidget parent=None, QString caption=QString(), QString directory=QString(), QString filter=QString(), QString encoding=QString()): argument 3 has unexpected type 'QPyNullVariant'

  Versione Python: 2.7.10 (default, Oct 23 2015, 18:05:06) [GCC 4.2.1 Compatible Apple LLVM 7.0.0 (clang-700.0.59.5)]
  Versione di QGIS: 2.12.1-Lyon Lyon, exported

I cannot enter the path manually. It does not allow me to click in the blank field.
Can I insert a code via the terminal?


texts:
i've come across an issue with the generated feature ID from GeoServer. The generated ID is a concatenation of the layer id and the feature primary key, rendering something like 'assets-1234566' say.

Is there a way to control this format of the ID?

I'm using vector tiles rendered with the Mapbox js SDK and some of the SDK functionality requires the feature ID, but it can only be an integer. The string returned from GeoServer is not accepted by Mapbox so certain functions do not work as a result.

Anyone struck this before? Any ideas on a work around?


texts:
Is there a way to update metadata of an existing geotiff file? I'm particularly interested in setting the nodata value. The following does not appear to work

import rasterio
from pprint import pprint

path = '/path/to/file.tif'

with rasterio.open(path, 'r+') as src:
    src.meta['nodata'] = -32767
    pprint(src.meta)


texts: How to create vector roads center lines from raster?

texts:
I need to know how to create roads vector center lines from raster using ArcGIS.
I have tried to use Iso Cluster Unsupervised classification tool to classify raster and set it to 30 classes. After that Reclassify tool to filter the roads cells from other values. After that I have used "Raster to polyline" tool and checked the Simplify box, but the results are a lot of lines looking like a polygons (screenshot). I need it to be just one road center line.
This is a classic topographic map and I have a lot of roads to process so this can`t be a semi-manual process, I need it automatically as much as that is possible.

Does anyone know better way to solve my problem?


texts: Adding Yandex Satellite map in OpeLayers

texts:
How can I add Yandex satellite map in OpenLayers?

new ol.layer.Tile({
  title: "Yandex Uydu",
  baseLayer: true,
  visible: true,
  source: new ol.source.XYZ({
    url: 'http://vec0{1-4}.maps.yandex.net/tiles?lang=tr_TR&l=s&x={x}&y={y}&z={z}',
    maxZoom: 19,
    attributions: [new ol.Attribution({
    html: '© Yandex'
    })],
     projection: 'EPSG:3395',
     tileGrid: ol.tilegrid.createXYZ({
       extent: [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]
     })
   })
})


texts: Mobile Template - ArcGIS API for JavaScript

texts:
I added dynamic map service to the mobile template from ESRI.
However, I see a couple of issues:

map does not cover full screen (when vertical).
background shows black on mobile browser (chrome).

Here is a site.


texts: How to filter and download OSM data by attribute tag?

texts:
Is there any online service that allows to download only one thematic layer of the OSM data for a specific geographic region? For example all objects related to vegetation (trees, grassland,...), but not streets, railways etc. There must be an easy and straightforward way to do so (I hope?).
I'd prefer to use the data directly in QGIS.

I am aware of Geofabrik which provides some data packages. But I am interested if there is any way to specify the data I want before actually downloading it?


texts: Lining attributes up with features when creating shapefile using PyShp?

texts:
I used the Python Module Pyshp to select particular census blocks from the Census's shapefiles and package these into a single shapefile. However when I add this data into ArcMap and right click to identify blocks, I get attributes for multiple blocks nowhere near the intended block. At first I assumed I had written my script in a way that misaligned the records and shapes objects in pyshp, but when I read this shapefile back into python, everything appears to be aligned correctly (i.e. I can look up the record I want and find the correct shape at the same index and vice versa).

Does anyone have a suggestion as to the next place I can look to try and solve my problem?

Here is the class I built to do this. The order of operations is to call load to get a particular zip file from the census, use add_blocks to collect any that match my list and add them to the writer object, and then write_shapefile to save whatever I've collected.

class CensusConnection:

    def __init__(self):
        self.file_name = ''
        self.ftp_path = 'ftp://ftp2.census.gov/geo/tiger/TIGER2010/TABBLOCK/2010/'
        self.reader = None
        self.writer = shapefile.Writer()

    def load(self,state_fips, county_fips = ''):
        '''Use me to load all blocks into a folder in pwd that match state FIPS, or state FIPS plus county FIPS.'''
        #import data and set up reader object
        self.file_name = 'tl_2010_' + state_fips + county_fips + '_tabblock10'
        tempdir = tempfile.mkdtemp()
        if not os.path.exists('tl_2010_' + state_fips + county_fips + '_tabblock10'):
            print 'Retrieving ' + self.ftp_path + self.file_name + '.zip...'
            sock = urllib2.urlopen(self.ftp_path + self.file_name + '.zip')
            memfile = io.BytesIO(sock.read())
            ZipFile(memfile, 'a').extractall(tempdir + "/" + self.file_name)
        self.reader = shapefile.Reader(tempdir + "/" + self.file_name + "/" + self.file_name)
        self._dump()
        #set up writer object to match reader
        self.shapes = self.reader.shapes()
        self.fields = self.reader.fields
        self.records = self.reader.records()
        self.writer.fields = self.fields

    def _dump(self):
        if (os.path.exists(self.file_name)):
            shutil.rmtree(self.file_name)

    def add_blocks(self, block_keys):
        '''Take a block key, find that block and add shapefiles to the writer object'''
        '''Block Key = StateFIPS + CountyFIPS + TractID + GroupID + BlockID (15 digit)'''
        fields = self.reader.fields
        records = self.reader.records()
        for index in range(len(records)):
            if(records[index][4] in set(block_keys)):
                self.writer.shapes().append(self.reader.shapes()[index])
                self.writer.records.append(self.reader.records()[index])

    def write_shapefile(self,dir_name):
        if(len(self.writer.records)>0):
            if(os.path.exists(dir_name)):
                shutil.rmtree(dir_name)
            os.mkdir(dir_name)
            self.writer.save(dir_name+"/"+dir_name)


texts: Removing overshots (dangles) using FME or FME/MRF Transformers

texts: I have a network of many connected lines where there is overshoot topological errors. Is there any way to remove them using FME Desktop?


texts: Identifying which shared boundary is longest using QGIS

texts:
This is an evolution from How to specify dependable layer names in QGIS Processing Modeller but should be considered completely separately.
I am creating a model in QGIS (3.10.3) Processing Modeller (Win 10[64]).
Final Objective
To create Voronoi polygons from Enumeration District (ED) centroids that are clipped to a higher, polygon, geometry.
Hurdle
The clipped Voronoi geography contains instances that either do (Voronoi-with-centroid) or do not ('merger-candidates') contain an ED centroid.
Chosen method to get over the hurdle
I create a unique FragmentID for every ED Voronoi instance (whether it has a centroid or not).
Then identify the neighbours of every ED Voronoi instance, and the shared boundary lengths as per:
https://spatialthoughts.com/2020/04/08/calculating-shared-border-lengths/ (by SpatialThoughts' Ujaval).
I now have a line (MultiLineString) dataset.
The next steps are (for each 'merger candidate') to identify their neighbouring 'Voronoi-with-centroid' with which it shares the longest boundary - and is in a different district.
The fields of interest are:
1: ZONECODE e41110005001  - the code of the ED Voronoi of which the 'merger candidate' fragment was a part (before clipping).
2: FragmentID e41110005001_67  - the code of the 'merger candidate' (after clipping).
3: NEIGHBOURS e41110109003_8,e41110005003_14,e41110005001_15,e41110109003_14,e41110005003_66,e41110108001_68,e41110110002_71
The codes of the Voronois ED fragments with which this ED 'merger candidate' fragment is a neighbour.
4: LENGTHS
,,2340,148,83,2289,381
The lengths of those boundaries
The BOLD parts of the CODES are the higher-level geography, the districts, to which the EDs belong. The Italic parts of the code are the EDs.
So, how would I identify that in this case, the neighbour with which e41110005001_67 shares the longest boundary, and with which it should be merged, is e41110008001_68 (2,289m)?
Yellow area e41110005001_67 need to be allocated to e41110008001_68


texts:
I run large site focused on history. Our current project expects the use of maps (for simplicity, we use current maps, not historical) with quite many points (custom SVG markers), lines, polygons on it with date (and if possible, time) parameter. I need a tool to create such items on existing map layer, to be able to place this map with features in form of slippy map on my website (without editing, just viewing), and to link clicks on the objects to some function that will load additional data (by some value field taken from the geodatabase).

So far, I was experimenting with QGIS desktop, where I was able to achieve
- use background map layer (via OpenLayers map plugin, Google, OpenStreepMap data)
- put on SVGmarkers, lines and polygons, save locally in ?shapefile?
- add value fields and date to marker table

So let's say, I am able to CREATE the features on the map.

But, I don't like the look of the background map layer, as Google is not showing the map 100% correctly in areas where I need it, on the other hand OpenStreetMap is great but it is too detailed for me (main roads, cities with names, water/rivers, country borders would be sufficient).

I'd love to have height shading added to the map, to create look similar to Google Physical map.

While reading a tons of information on the internet, I came to following requirements:

I need a geo app server with database - doable with GeoServer (Java) and PostGIS
I need to grab data from geofabrik (I plan to use it for the whole world with most of our focus will be Europe, Asia, Africa and North America) and load it into PostGIS (using tool like osm2pgsql)
I need to render the tiles from this data with tool like Mapnik
and I need to publish my layer from QGIS to GeoServer, then use some JavaScript library (Leaflet) to show it on my page.

Am I right?

And how do I pick only those information from the OpenStreetMap that I need on my map? I expect it to look similar to this

http://www.geofabrik.de/maps/rendering.html
Geofabrik TOPO STYLE

Footnote: I do not need any special GIS measurement tools, geocoding, path finding, no no. Just map and items to be shown on the map, in various layers (able to turn them on and off in browser), browsable by date, with possible thousands of objects all around the globe. I cannot predict the usage, hardware is no limit for me, so I prefer to create the full stack on our server and run it from there, not to pay for some ready service (that might or might not fully fit into our needs).


texts: MapInfo Professional: Can't modify table structure

texts:
I am new to MapInfo and I need to combine several tables together, by copying and pasting.  I need for all the tables to have the same table structure, however when I try to modify these tables (table>maintenance>table structure>view/modify table) it won't allow me to change the order of the fields.  When i try to, i get the following error: "cannot perform this operation on table (table name)." Why won't it let me change the data structure? any suggestions?

The picture is of the dialogue window i get when i attempt to modify table structure.  As you can see, it doesn't look like the regular window.  But i don't know why or how to change it.
When i move my fields to the order I need them in, and click OK, is when i get the error message.

THanks!


texts:
When I load a .vrt file, the image comes in with a black and blue tint usually or is some other weird color.  If I zoom in at all QGIS crashes and creates a minidump in my program folder. Also, when I try to Save-As the .vrt file (as a .Tiff and .vrt) the program entirely crashes and doesn't respond.  These are multiband rasters.

My ultimate goal is to turn the .vrt into a .tiff file.

I have all the images I used to create the .vrt file in case there is another route to accomplish this.


texts: limit photon geocoding results to a specific place

texts:
I would like to geolocate addresses from a particular city.

I have just recently used photon R package as a geolocation solution with certain degree of success. Unfortunately, some results are placed very far away from the place I am trying to geolocate, even if I use location bias:

df = geocode(address, limit = 1,key = "highway", locbias = c(-0.8787, 41.6486))

As far as I know Photon can accept keys and values, but don't know how to use it, further than what I did (If I am not wrong it will search only OSM's objects tagged as highway=*). Is there a way to tell phonon to search only OSM's objects tagged as highway=* AND are located within a city (maybe using OSM's is_in=*)?

Reproducible code:

# Only the last address is located in Zaragoza (OK). The others are not (KO)
address = c("Pablo Neruda, Grancasa, Zaragoza", "GAVIN, 14, La Magdalena, Zaragoza", "La Magdalena, Zaragoza", "calle Ciprés, 8, Alfonso, Zaragoza", "calle asalto, 52, La Magdalena, Zaragoza")

df = geocode(address, limit = 1,key = "highway", locbias = c(-0.8787, 41.6486))

Note: from what I see, some wrong results are due to the fact that are not proper streets but neighbourhood names, hence my key = highway is messing things up. Still, some others are streets that do exist in Zaragoza, although with some minor spelling differences -hence that's why Mexico scored better).


texts: Setting the column-order of the attributetable in ArcMap programmatically?

texts:
I'm loading a feature layer to ArcMap and want to prepare the columns shown in the attribute-table.
Aliasname and visibility I set with ILayerFields-Interface of the FeatureLayer by passing it to the IFieldInfo. But I haven't found a way to set the order of the columns.

Found out reading the column-order is possible with FindFields:

public virtual object FindFields { get; }
    Member of ESRI.ArcGIS.Carto.FeatureLayerClass

It gives me the column-names in choosed order. But it has only a get and not a set.

Does anybody has an idea?


texts: How to create a point between two other points in QGIS?

texts: We are capturing sample points from survey maps. In between two points we have to introduce another point (as in fig), like - if the sample points are 42 and 42.5 we have to introduce another point with value 42.5. Is it possible to automate this process in QGIS using any plugin ?


texts: Locating property corners using plat map description?

texts:
I am interested in a piece of property for which I have a plat map description, giving what seem to be latitudes/longitudes of lines along the boundaries (not sure about this). The objective is to locate the corners of the rectangular parcel.

The data given seems to define some lines (not points), and is in reference to portions of sections. It seems that I need to know the lat/long of the section corners to do this, but I am unable to find them.

How can I find the corners of the section?

Here is part of the description:

  ... sw 1/4 of ne 1/4 , sec 3, twshp 25 North, Range 2 East, W.M., in
  Kitsap County, WA, described as follows: Beginning at a point North 89
  deg. 49' 02" West 71 ft. FROM THE SE CORNER OF SAID SW 1/4 OF THE NE
  1/4; thence North 7 deg. 7'04" East 440.42 ft. ...

I want to plug lat/long into my GPS and walk to the corners.


texts: Creating Workspace, Data Store and Layer via GeoServer REST API?

texts:
I'm trying to automate the process of setting up a Workspace, PostGIS Data Store and Layer in GeoServer using the GeoServer REST API.

Using the excellent answer in Create a Layer in GeoServer using REST I'm able to do this, with one sticking point. When I create the workspace using this code...

import requests, json
headers = {'Content-Type': 'application/json'}
auth = ('admin', 'geoserver')
url = "http://<my_geoserver>:8080/geoserver/rest/workspaces"
data = {"workspace": {"name": "blah"}}
r = requests.post(url, headers=headers, auth=auth, data=json.dumps(data))

... the Services are all unchecked on the newly created Workspace:

According to How to set up a Workspace setting in Geoserver using Curl this is a missing feature in the GeoServer REST API.

Given that the Services are unchecked, I'm unable to publish a layer via the GeoServer REST API:

url = 'http://<my_geoserver>:8080/geoserver/rest/workspaces/<workspaceName>/datastores/<dataStoreName>/featuretypes/'
data = {"featureType": {"name": "<layerName>","srs": "EPSG:4326","enabled": "true","store": {"@class": "dataStore","name": "<workspaceName>:<dataStoreName>"}}}
r = requests.post(url, headers=headers, auth=auth, data=json.dumps(data))

as the API responds with :Schema 'http://<workspaceName>:<layerName>' does not exist.

If I manually select the 4 Services checkboxes on the newly created Workspace, then re-run the code to create a layer, the Layer is created without error.

Is it possible to completely automate the process of creating a Workspace and Datastore, then creating Layers within that Datastore and Workspace, via the GeoServer REST API?


texts:
What's a proper data type to save coordinates in QGIS when using the 'Move label' tool?

real
double
numeric()
...

EDIT:

The CRS of the project is MGI / Austria GK East (31256).


texts:
I have a layer that describe rainfall for 10 years. Now I want to interpolate the area in ArcGIS Desktop (ArcMap) 10.7. As you know I should define some point on my area to extract the rainfall from the Layer.

How should I define these points i.e. in the center of pixels or more?

This is my Layer and a closed Border:


texts: Finding lat long in ArcGIS similar to as in Google Earth

texts:
I want to know if there is some possibility to find lat long position in ArcGIS 10.0 as we can find on Google Earth. I have an assignment and doing this can save me a lot of time. I have already used the "Find" option in ArcGIS 10.0 which is limited to what I want.

In the question above use of "Google Earth" have been taken as an example only. I would like to find the lat long position in the platform of ArcGIS as I have been trying to find with "Find" tool, which did not help.


texts: How to create shadows of buildings with r.sunmask?

texts: I need to create image-raster with shadow of relief and shadow of buildings in one. I know how to create shadow of relief but how i can create shadow of buildings with r.sunmask? I dont know how from 2D polygon layer of buildings create shadow of complete buildings.


texts: ArcGIS REST layer doesn't show up in OpenLayers

texts:
I try to show a WMS layer with OpenLayers 4.6.5 but the tiles don't show up.

I use the following code to load the WMS layer as ArcGIS REST service:

new ol.layer.Tile({
  source: new ol.source.TileArcGISRest({
    url: "https://geoportal.amb.cat/arcgis/rest/services/Cartografia/MapaTopograficAMB1M_ETRS89/MapServer",
    projection: 'EPSG:25831',
    params: {'LAYERS': 'Nivell 1k'}
  })
})

The browser network tab shows that tiles get loaded, but they are all transparent.

Here is the jsfiddle.


texts: Is it possible to import SAGA colour styles in QGIS?

texts:
iam still using QGIS and SAGA independently in my productive use (despite the excellent SEXTANTE toolbox, which is sometimes a little buggy).
SAGA has some very nice default styles and i am wondering if it is possible to import SAGA color palettes (extension .pal) into QGIS?


texts: Why Lat and Long values through attribute table in ArcMap are not accurate?

texts:
I calculated the lat, long in metres through the attribute table, and the values are very course (i.e., varying in +/-50m only).

But, when I move the mouse cursor over the points, the lat long values are much finer (mm)

Example: It's varying in around 22 metres with x coordinate when I moved from 1 point to its immediate neighbor point, but it's not reflected in the attribute table.

Refer the screenshots in the attachments.

How can I tackle this?


texts: Adding attributes to CSV-file using PyQGIS

texts:
I am new to PyQGIS and I want to add one more field in my existed table. Here is my code:
pointLayer = QgsVectorLayer(uri, 'my_layer', 'delimitedtext')
pointLayer.dataProvider().addAttributes([QgsField("new1", QVariant.Double)])
pointLayer.updateFields()

However, there's no difference if I print the names of the field before and after adding attributes respectively.

I don't know how to fix the problem and the code does not have any error.


texts: How do I use GDAL to convert from an arbitrary projection/datum into Equirectangular/WGS84?

texts:
I've been trying to use GDAL to convert from an arbitrary projection/datum into Equirectangular (aka Plate Caree, aka Equidistant Cylindrical, aka Geographic, aka Lat/Long) with WGS84 datum.  All of my input files have an associated PRJ file that I've been using with GDAL with this statement:

gdalwarp -s_srs Input.prj -t_srs Output.prj Input.xyz Output.tif

So far, I've had mixed results.  I've been able to successfully transform a Lambert Equal Area projection in WGS84 to Geographic/WGS84, but performing a transform on LEA in RT90 projection causes an error of ~1km.

Stranger still, I tried to transform two Sinusoidal projections in RT90 and in WGS84, both of which resulted in an error of 16km!  I suspect that these errors are largely due to some sort of ellipsoid/sphere discrepancy.

I'm very confused as to why these errors are so random, but more concerned that they are happening at all.  Shouldn't gdalwarp take any ellipsoid information about the source and target from the PRJ files and transform between them?

What can I do to correct my issue?


texts: How to fill attribute field of one feature based on calculation on attribute field of another feature

texts:
I would like to fill an Attribute field of a geodateabase table based on some calculations such as 'sum' over specific field of another Feature class. Is there an Approach to formulate it? while in the "field calculator" we are able just to introduce a Formular inside that specfic Feature.

I am thinking about coding for this reason but to be honest as I am not pro in coding, I am sure what is the right way to think about it.

I would be thankful to have your assistant with my python code.

   import arcpy

   from arcpy import env

   # Set the current workspace
   env.workspace = "K:\Lines\M966KAZ\55_SdiGis\45_Gdb\File\Test_Corridor_relationship.gdb"

  table = "link/to/Corridor_stat.dbf"
  fc = "link/to/Corridor_test_1.shp"

  # Create a search cursor
  rowsFc = arcpy.SearchCursor(fc)

  # Create an update cursor
  rowsTable = arcpy.UpdateCursor(table)

  for rowT in rowsTable:
     if fc.Corridor_S <= 1:
       valueTable = sum (fc.SHAPE_Leng)
  elif (fc.Corridor_S > 1 and fc.Corridor_S <= 60):
       valueTable = sum (fc.SHAPE_Leng)
  elif (fc.Corridor_S > 60 and fc.Corridor_S <= 100):
      valueTable = sum (fc.SHAPE_Leng)


texts: Main Menu and Catalog Tree View in ArcCatalog are not visible

texts: In ArcCatalog (the program, not the window in ArcMap), the Main Menu (File, Edit, View... Windows, Help) at the top and the Catalog Tree View have disappeared and I can't find a way to reconnect them. How do I make them visible again?


texts: Converting Polygon into single or all Points using QGIS?

texts:
I want to convert polygons into points using QGIS.

I know how it's done in ArcGIS Desktop.


texts: Setting region grow threshold by layer value difference using eCognition?

texts:
I'm trying to use the "grow region"-algorithm in eCognition Developer 9 on an yet almost unclassified 1x1 (=pixel resolution) chessboard segmentation. I classified nothing but my "seed points" by a thematic layer yet.

I have a single layer that I want to use for growing ("infinite cycles", as much as possible) as long as the candidate neighbors value in this layer is within a specific distance.

f.e.: let's say i don't want to add a neighbor if the value difference exceeds +-2.

How could i do this within a ruleset?


texts: Spatial reference does not match data frame for Bing Maps

texts:
In ArcObjects I am loading Bing Maps as a base layer. It's coordinate system is esriSRProjCS3Type.esriSRProjCS_WGS1984WebMercatorMajorAuxSphere. On top of the base layer, I am trying to load points, lines and polygons via SqlGeometry objects where geometric coordinate system is WGS1984 or esriSRGeoCSType.esriSRGeoCS_WGS1984.

I am converting SqlGeometry to ESRI's IGeometry objects and set spatial reference to esriSRGeoCSType.esriSRGeoCS_WGS1984 specifically:

public static IGeometry Convert(SqlGeometry sqlGeometry) {
  int count;
  IGeometry esriGeometry;
  var factory = (IGeometryFactory3) new GeometryEnvironment();
  factory.CreateGeometryFromWkbVariant(sqlGeometry.STAsBinary().Value, out esriGeometry, out count);
  esriGeometry.SpatialReference = spatialReferenceForGeoCS_WGS1984;
  return esriGeometry;
}

And I am following ESRI's documentation to create a feature class and assign spatial reference to it (doc):

public IFeatureClass CreateFeatureClass(String featureClassName,
  IFeatureWorkspace featureWorkspace,
  ISpatialReference spatialReference)
{
  // Instantiate a feature class description to get the required fields.
  IFeatureClassDescription fcDescription = new FeatureClassDescriptionClass();
  IObjectClassDescription ocDescription = (IObjectClassDescription)fcDescription;
  IFields fields = ocDescription.RequiredFields;

  // Find the shape field in the required fields and modify its GeometryDef to
  // use point geometry and to set the spatial reference.
  int shapeFieldIndex = fields.FindField(fcDescription.ShapeFieldName);
  IField field = fields.get_Field(shapeFieldIndex);
  IGeometryDef geometryDef = field.GeometryDef;
  IGeometryDefEdit geometryDefEdit = (IGeometryDefEdit)geometryDef;
  geometryDefEdit.GeometryType_2 = esriGeometryType.esriGeometryPoint;
  geometryDefEdit.SpatialReference_2 = spatialReferenceForGeoCS_WGS1984;

  // Create the feature class.
  IFeatureClass featureClass = featureWorkspace.CreateFeatureClass(featureClassName, fields,
    ocDescription.InstanceCLSID, ocDescription.ClassExtensionCLSID, esriFeatureType.esriFTSimple,
    fcDescription.ShapeFieldName, "");
  return featureClass;
}

Still when I start editing session, I am getting "Spatial reference does not match data frame" error. Is there another artifact that needs spatial reference?


texts: Creating Lines from Multiple Points to one Reference Point

texts:
I'm looking to create a new line shapefile to connect two sets of point files using a common reference field.

E.g. in the below image i'd have 7 lines connecting from the blue shapefile to the orange shapefile. I have a few hundred clusters of features like this so need a fairly automated way.


texts: Looping through AGOL services with arcrest

texts:
Using this example I am trying to loop through my Arcgis Online account using ArcREST in order to download feature layers locally.  However I have not been able to drill down into each folder to acquire each service name with .content in order to run the .exportItem.  I can get the folder names but I can't get access to the feature layers in the folders.

import os
import arcpy as env
import arcrest
import arcpy

arcpy.env.workspace = r"C:\xxx\temp"
uname = "xxx"
pword = "xxx"

site = "http://xxx.xxx.arcgis.com"

sh = arcrest.AGOLTokenSecurityHandler(uname,pword,site)
admin = arcrest.manageorg.administration.Administration(securityHandler=sh)

content = admin.content
user = content.users.user(uname)

for folder in user.folders:
    print folder
    #This line below was added after original post
    for sub_folder in folder:
       print sub_folder

Edit: added a for loop within the original loop to attempt to access the contents of the folders.  The variable sub_folder prints only the contents of the root folder


texts: How to view a postgis raster query in QGIS

texts:
I have successfully loaded a large number of .asc images (>2800 200x200 tiles) into a Postgres table with PostGIS 2.1 support. The new table appears in QGIS's DB Manager and I can add the whole image as a single layer.

Is there a way, for example, to create an SQL script on raster tables and display the result in QGIS as a raster layer as you might with a vector layer? Or simply to filter the table to specific tiles.

Within QGIS, the DB Manager allows me to load the whole table but doesn't seem to support filtering of the table, and neither does right clicking on the layer once its loaded.

This is my first foray into PostGIS and rasters so the answer may be obvious so please have patience with me.


texts: Using Arcpy.Mapping to add raster layers?

texts:
I am trying to add raster layers in specific data frames @current mxd. The script works in the ArcGIS 10.3 Python console but when I created a custom toolbox and add the script in the custom toolbox the script tool neither works nor gives an error.

Here is the script:

import arcpy,os
from arcpy import env

mxd=arcpy.mapping.MapDocument('Current')
workspace=arcpy.GetParameterAsText(0)
dflist=arcpy.mapping.ListDataFrames(mxd)
listraster=arcpy.ListRasters()
for raster in listraster:
    lyr=arcpy.mapping.Layer(raster)
    if str(lyr)=='om_in_percent':
        arcpy.mapping.AddLayer(dflist[0],lyr)
    elif str(lyr)=='cec_status':
        arcpy.mapping.AddLayer(dflist[1],lyr)
    elif str(lyr)=='n_inpercent':
        arcpy.mapping.AddLayer(dflist[2],lyr)
    elif str(lyr)=='si_status':
        arcpy.mapping.AddLayer(dflist[3],lyr)
    else:
        print(arcpy.GetMessage(1))

Any Help?


texts:
I combined two shapefiles and noticed that the polygons didn't overlay perfectly.
How can I resolve these intersections to snap to their closest neighbours?

This is a closer look at the intersections:

I brought these files together using mysql, however I've been using ogr2ogr and python also.


texts: Using time manager for shapefile with different timestamps

texts: I have a network which contains speed information for each segment at different timestamps. I would like to create an animation of speed changes during the day. There is multiple columns which contain speed for a segment at a given time. How do I work around the fact that time manager plugin only reads time value from one column ?


texts: Convert batch file using GDAL_Translate.

texts:
how to convert batch file from .tif to .bil format using gdal_translate in windows CMD command line, file in the directory.

Directory path : D:\GDAL-TEST
001_mean.tif
002_mean.tif
003_mean.tif
004_mean.tif
So on....

I have one code but by this code we can convert single file at a time but i want to loop through all raster which have .tif extension in the  directory and convert it to .bil format. Please want a solution

My code. I a running this code on GDAL Command window on windows 7 platform.

gdal_translate -of "ENVI" D:\GDAL-TEST\001_mean.tif D:\GDAL-TEST\001_mean.bil


texts:
I've been trying to create a 400m buffer in a polygon layer in QGIS 2.18 and it seems like I may be doing something wrong compared to how the tool used to work.

Does the tool still use the units of the CRS?

When I try to create a 400m buffer or even a 10m buffer I end up just getting a large buffer that definitely isn't either one of those distances.   I'm going to Vector->Geoprocessing Tools->Fixed Distance Buffer and setting things as follows in the screen shot.  GeoJSON is in EPSG2768.

Seems like the buffer tool wasn't using the defined CRS units for my layers.  Checked the output buffer layer and it was using WGS84 instead of EPSG2768.  CRS for the project and all the layers was set to EPSG2768.  Things didn't work until I changed the CRS for the buffer layer and then re-ran the buffer tool.  I set up a couple virtual machines to test things out with mixed results.  Fresh Xunbuntu VM with QGIS installed and the buffer tool worked right off the bat, but a fresh KDE Neon VM with QGIS had the same issue I initially posted about.  Things are working now though.


texts:
I want to add WFS layer from MapServer in an OpenLayers 4.0.1 client.

  var mapserver_wfs_source = new ol.source.Vector({
    format: new ol.format.WFS(),
    loader: function(extent, resolution, projection){
      var url = 'http://localhost:8081/cgi-bin/mapserv.exe?map=/ms4w/apps/tutorial/htdocs/example2_wfs.map&SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=hydro';
      $.ajax({
        url
      }).done(loadFeatures);
    },
    strategy: ol.loadingstrategy.bbox,
    projection: 'EPSG:4326'
  });

  var loadFeatures = function(response) {
    //console.log(response);
    mapserver_wfs_source.addFeatures(new ol.format.WFS().readFeatures(response));
    console.log(mapserver_wfs_source.getFeatures());
  };

  var mapserver_wfs = new ol.layer.Vector({
    source: mapserver_wfs_source,
    style: new ol.style.Style({
      stroke: new ol.style.Stroke({
        color: 'rgba(0, 255, 255, 1.0)',
        width: 2
      }),
      fill: new ol.style.Fill({
        color: 'rgba(0, 0, 255, 1.0)'
      })
    })
  });

  var projection = new ol.proj.Projection({
    code: 'EPSG:4326',
    units: 'degrees'
  });

  var map = new ol.Map({
    layers: [osm, mapserver_wfs],
    target: 'map',
    view: new ol.View({
      center: [-90, 45],
      zoom: 4,
      projection: projection
    })
  });

If I just follow the link 'http://localhost:8081/cgi-bin/mapserv.exe?map=/ms4w/apps/tutorial/htdocs/example2_wfs.map&SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=hydro' I can see the features, and "mapserver_wfs_source.getFeatures()" returns the features.

Any idea on what is going on here?


texts: Converting shapefile to WKB after v.clean

texts:
I've been trying to get to through this by reading manuals, but I can't seem to understand what I'm doing.

I have source data with overlaps that exist in a table with associated WKBs (SRID 4326).  I eliminated these overlaps using the v.clean tool in the processing toolbox panel.  The output that I receive is a shapefile and at first glance seems to overlap perfectly (the overlaps are resolved when run against the topology checker).

When I run Import that shapefile into PG as a table via Shapefile Loader 2.1 plug-in and compare the visualization it is a very different scale and location as shown in the image below.

The original image is pictured on the right.  Perhaps I'm approaching this problem the wrong way, but it seems simple enough to grasp what I'm trying to do. Ideally I maintain the same dimensions the above image but modify the underlying source data.


texts: QGIS distribution with automatic installation of plugins

texts:
in our company we use QGIS 3.4.8LTR. We install qgis on every PC via distribution software. We manage the user settings with the qgis_global_settings.ini. So we can ensure that every user has the same settings.
Now we are searching for a possibility to manage the installtion of plugins. With the update on a new QGIS version we want, that every user has a pre-selection of plugins automatically installed.
How can we manage that?


texts: Publish graphic icon into geoserver from QGIS

texts:
I want to publish graphic icon (point) into geoserver by using maki style. Below is the coding,

<UserStyle>
  <se:Name>streetlight_u4</se:Name>
  <se:FeatureTypeStyle>
    <se:Rule>
      <se:Name>Single symbol</se:Name>
      <se:PointSymbolizer>
        <se:Graphic>
          <se:ExternalGraphic>
            <se:OnlineResource xlink:type="simple" xlink:href="file:///maki/park-alt1-15.svg"/>
            <se:Format>image/svg+xml</se:Format>
          </se:ExternalGraphic>

But, the result turns like below:-

Which is the actual graphic point is like below :-

Can anyone help?


texts:
I'm creating an Esri Story Map and trying to set a fixed extent for my maps.  For example, if I'm trying to illustrate fire stations in the City of Denver, CO I don't want my users to be able to zoom/pan half way across the globe on accident.

I understand that I can set an initial extent - http://www.esri.com/esri-news/arcuser/fall-2013/control-the-map-extent-in-arcgis-online but this doesn't prevent the user from zooming/panning away from my data.  I tried to customize my URL according to this source https://doc.arcgis.com/en/arcgis-online/reference/use-url-parameters.htm but haven't been successful.


texts: Obtain voronoi polygons which extend till the boundary of another polygon

texts:
I have certain hubs (coloured yellow) in a district (coloured light blue). I need to find polygons indicating the closest hub to every point in the district.

The voronoi polygons that I obtained using QGIS are coloured brown in the image. However, I need the voronoi polygons to fill the outer boundary as I need the polygons to fill the entire district.

Adding a buffer region while making the voronoi polygons and then using the intersect function to crop the polygons to the district boundary gives me wrong polygons as shown by the following image:

How do I generate the voronoi polygons that I require?


texts: Maximum likelihood image classification in ArcMap

texts: I am running MLC in ArcMap but the result is returning 1 class instead of all 3 classes that the training samples have. What could be the cause of that? is my sample too little or what could cause it?


texts: Read SDE Feature Attributes without Arcpy

texts:
I've been building some software that reads/writes spatial and non-spatial data out of our SQL SDEs dynamically when called. As I've pointed out before Time Delay in Importing ArcPy , ArcPy takes forever to import. This has become an issue with my users as the querying and writing tools are too slow for their day to day use.

I am reading and writing data to tables and spatial features within the SDEs using ArcPy cursors (search, update). I've found that cursors work best for the reading and writing of data within SDEs, but I am hoping there is a quicker way without ArcPy.


texts: Create grid points and compute distance within regions

texts:
The spatial analysis problem is as follows. Suppose you want to compute to average distance a consumer within each postal code has to travel to get the nearest bank. I have coordinates of bank locations (using the sf package, the geometry is point). And I have a shape file with postal codes (this is given as a multipolygon geometry type).

At this stage I want to create a grid with very fine points and overlay it on top of the postal code shape file. Using these points, I need to find the mean distance each point as to travel to get to the nearest bank location. The idea then is the to compute the average by zip code.

I have it conceptually worked out, but I'm stuck on the part where I have to make a fine grid. I'm not even sure how to approach this given I have two difference data sets with different geometries. I'm using the sf package and the tidyverse.


texts: external measurement connection to qgis experimental

texts: Recently I used QGIS experimental (for Android). I need a direct connection from a mobile or tablet measurement device to QGIS experimental. I added external resources field for feature photo which is expected to be taken directly on the field by a mobile camera and located via GPS. Is this possible in QGIS experimental (not Qfield)?


texts: GeoServer - compute from native bounds does not work correctly

texts:
I am currently working on the topic GeoServer, so I am very new.
For some layers it is not possible to calculate the BBox via the button "compute from native bounds".
What can be the problem?
What is the best way to troubleshoot?


texts: GDAL's ReadAsArray() incredibly slow when reading an in-memory raster?

texts:
So I have a script set up in a conda environment on macOS Catalina running python 3.7, I'll post a paraphrased version of it:
from osgeo import gdal
import numpy as np

grib_url = '/vsicurl/some_grib_url'
grib_ds = gdal.Open(grib_url)

translated_grib = gdal.Translate('', grib, **translate_kwargs)

reprojected_grib = gdal.Warp('', translated_grib, **warp_kwargs)

band1_arr = reprojected_grib.GetRasterBand(1).ReadAsArray()

Essentially, I'm getting a grib file from a URL, reading it in memory, cropping and retrojecting it in memory, pulling out a specific band, then converting it to a numpy array. However, the script will not finish in the terminal. It's been hung up for a couple hours without throwing any errors.
If I convert the band to an array before translating and retrojecting the entire grib, it will convert it fine (within about a minute).
I've tried setting environments with python 2.6, 2.7, 3.6, 3.7, and 3.8 and the problem persists. I've also tried various combinations of gdal/numpy versions as I've seen that's a common problem. I've also tried iterations of installing one before the other.
I don't think there's necessarily an error, but rather gdal is just really slow opening a VRT in memory? Has anyone had a problem similar to this or know a way to fix it?


texts: Reducing ee.Image to a list and then exporting as CSV file

texts:
I am struggling with transforming a single image to a list and then exporting it to my drive.

I am using the ee.Reducer.toList() but the output looks like this:

Insted of this:

When I print the featureCollection that is to be written to the file - it seems to contain all the data needed.

Here is my code:

var region = ee.Geometry.Polygon([
  [[-5, 40], [-5, 60], [-5, 60]]
]);

//Dates of Interest
var start = ee.Date("2018-04-01");
var finish = ee.Date("2018-05-01");

var sen2 = ee.ImageCollection('COPERNICUS/S2_SR')
  .filterDate(start,finish)
  .filterBounds(region)
  .map(function(image){
    return image
    .clip(region)
    .normalizedDifference(['B8','B4'])
    .rename('NDVI')
  })

print(sen2)

var sngl = ee.Image(sen2.first())

var list = sngl.reduceRegion({
  reducer: ee.Reducer.toList(),
  geometry: region,
  maxPixels: 1e8,
  scale: 10
});

// Make a feature without geometry and set the properties to the dictionary of means.
var feature = ee.Feature(null, list);

// Wrap the Feature in a FeatureCollection for export.
var featureCollection = ee.FeatureCollection([feature]);

print(featureCollection)

// Export the FeatureCollection.
Export.table.toDrive({
  collection: featureCollection,
  fileNamePrefix: 'field_test',
  folder: 'field_test',
  fileFormat: 'CSV'
});

Any ideas?


texts: Geoserver SLD: How to get CENTROID of result of INTERSECTION?

texts:
In my GeoServer SLD I need to place a symbol (particularly a piechart) exactly on Centroid of the result of Intersection of 2 geometry:

geom field from my PostGIS table AND
environment variable WMS_BBOX

My aim is to place my Piechart on centroid of shown part of a polygon (an not on centroid of the whole polygon itself).

Ian Turton commented on my previous question with the following Geometry element in SLD, and to use it in <PointSymbolizer>:

<Geometry>
  <ogc:Function name="intersection">
    <ogc:PropertyName>the_geom</ogc:PropertyName>
    <ogc:Function name="env">
      <ogc:Literal>wms_bbox</ogc:Literal>
    </ogc:Function>
   </ogc:Function>
</Geometry>

The problem is, with the above Geometry (of intersection between geom field and WMS_BBOX) the piechart is shown in every geometry element of my multipolygon table. And this is surely not desirable.

This effect is similar to forcing a <PointSymbolizer> Rule (below) to a multipolygon table.

<Rule>
<PointSymbolizer>
    <Graphic>
        <Mark>
            <WellKnownName>square</WellKnownName>
            <Fill><CssParameter name="fill">#FF0000</CssParameter>
            </Fill>
        </Mark>
        <Size>6</Size>
    </Graphic>
</PointSymbolizer>
</Rule>

The above style enforcement gives the below layer (these islands belong to 1 row of a multipolygon table).

Again this is not desirable.

So then i tried the <ogc:Function name="centroid">. Thinking that this is equivalent to a PostGIS command of

ST_CENTROID(
    ST_INTERSECTION(
        geometry_field, --of my multipolygon table
        WMS_BBOX
        ))

... then I wrote the following in my SLD:

<Geometry>
    <ogc:Function name="centroid">
        <ogc:Function name="intersection">
            <ogc:PropertyName>geom</ogc:PropertyName>
            <ogc:Function name="env">
                <ogc:Literal>wms_bbox</ogc:Literal>
            </ogc:Function>
        </ogc:Function>
    </ogc:Function>
</Geometry>

However my WMS preview throws this error:

code:"internalError:Rendering process failed
Filter Function problem for function centroid argument #0 - expected type Geometry

This is as if Centroid argument expects a geometry type but Intersection argument gives non-geometry type. The GeoServer Geometric functions reference however, while clearly states the type of arguments input, does not states clearly the type of the output.

Question:

why does centroid argument fail on my intersection?
how to solve it?

My playground:

GeoServer 2.14.3 (with chart extension)
Apache Tomcat 8.5.32 Server
Windows 10 Pro 64-bit on Xeon E5-1660 3GHz, 32GB memory
PostgreSQL 9.6.11, PostGIS 2.5


texts: How to convert OSM relation XML data to shapefile?

texts:
I would like to transform the XML data here: http://wiki.openstreetmap.org/wiki/Berlin/Transportation into shapefiles or load into PostGIS.

I think the thing I am missing is how these XML files, which contains relations, can be used to get the data they represent.

I have tried opening the JOSM link, opening the GPX files in QGIS, and directing to the XML file to osm2shp.jar - none of these routes has given me what I need.  The GPX files can be used to generate a shapefile, but when I 'Open Attributes Table' in QGIS all of the column values are NULL for all points.

I would appreciate a step by step description of how to work with this data in PostGIS or as a shapefile.


texts: No legend showing at all in ArcGIS Server print task (MapExport)

texts:
I am using some stock print templates like the 2 in this sample. The problem is that the legend doesn't show up at all. I found a few forum questions showing problems with the legend but none of them involved having no legend at all. My code looks like this ...

    var legendLayer = new LegendLayer();
    legendLayer.layerId = "Parcels_and_Zoning";
    legendLayer.subLayerIds = [4, 10];

The LegendLayer is included in the require section "esri/tasks/LegendLayer" and using typeof(legendLayer), it is an object. None of the samples that I have found show where to include space for the legend so maybe that is where things are going wrong. This next bit of my code is almost directly from the sample that I linked above ...

        var layouts = [{
        name: "Letter ANSI A Landscape",
        label: "Landscape (PDF)",
        format: "pdf",
        options: {
          legendLayers: [legendLayer], // empty array means no legend
          scalebarUnit: "Miles",
          titleText: printTitle + ", Landscape PDF"
        }
      }, {
        name: "Letter ANSI A Portrait",
        label: "Portrait (Image)",
        format: "jpg", // ETC. ETC.

Note the comment saying "// empty array means no legend". I took the legendLayers option out entirely and STILL get no legend. Does anyone have any idea what I am doing wrong? I considered that I may need to add a legend in the html but I only need it for printing (I have not fully explored this). Do I need to create a legend in some way other than the code I've provided (the samples don't indicate this at all and all seem to avoid having a legend)? Is it possibly a problem with the service I have created ? The screenshot shows the result with no legend and sublayer 10 is the pink and amber colored polygon layer that I was hoping to get in the legend.


texts:
I'm working in ArcMap 10 and can move around in R, though am not an expert.

I have a shp file and an xlsx file.
The shp file is 20,221 rows and includes a field GID (grid ID) which is not unique among rows because multiple polygons named with a PID (park ID) can have the same GID. PID-GID pairs are unique.

The xlsx file 230,895 rows and also has a GID field which is not unique among rows because for each GID there are 60 years worth of data, each year represented by a row. Year-GID pairs are unique.

How can I combine the two files so that for each PID-GID combo in the shp file I get a row for each GID-year combo in the xlsx file?

I've been playing with joins and query tables and can't get it yet.


texts: How to know contents of column names in GRASS?

texts:
I am using GRASS to extract a specific area from a vector map. Here is the script that I am running:

v.extract PHL_map out=makati where="AREA_NAME = 'Makati"

The problem is, I don't know what are the values of the column names. I can only view the column names but I can't view the values under the column names


texts: Enable Handler Class Methods to Yield Instance Attribute in Pyosmium

texts:
This question is regarding the package pyosmium specifically. I was just wondering if the following functionality is possible, and if not how it could be implemented.

I want to stream/yield certain instance attributes instead of updating them in-memory.

Currently we can do the following:

class Handler(osmium.SimpleHandler):
    def __init__(self):
        osmium.SimpleHandler.__init__(self)
        self.edge_and_nodes = []
    def way(self, w):
        self.edge_and_nodes.append({'edge_id': w.id,
                                        'nodes': [w.nodes[i].ref for i in range(len(w.nodes))]})
h = Handler()
h.apply_file("test.osm.pbf")
print("Edges and their connected nodes: {}".format(h.edge_and_nodes))

However, when dealing with large regions this is not scalable.

I would like a way of yielding a dictionary object that includes WayIds and related NodeIds (as well as tags, etc) for every WayObject. Is this possible?

I am looking for something like this:

class StreamHandler(osmium.SimpleHandler):
    def __init__(self):
        osmium.SimpleHandler.__init__(self)
        self.edge_and_nodes = []
    def way(self, w):
        yield {'edge_id': w.id,
               'nodes': [w.nodes[i].ref for i in range(len(w.nodes))]}
h = StreamHandler()
h.apply_file("test.osm.pbf")
for row in h.way(w):
    print(row)

But I am not sure how to pass the w parameter (WayObject) since that seems to be dealt with internally using the apply_file() method (and I can't seem to find the source code for that method).


texts: Import PostGis tables into a GeoDataframe unexpected : 'AttributeError: 'NoneType' object has no attribute 'encode''

texts:
When I try to import two of my tables from postgis to geopandas I have an error. I guess the problem is due to my geometry columun, but I don't know why.

engine = create_engine('postgresql+psycopg2://xxx:xx@host/db')
    sql = 'Select * FROM table'
    GeoDF= gpd.GeoDataFrame.from_postgis(sql, engine)

It works great for most of my tables but two tables can't be imported, I have the message :  'AttributeError: 'NoneType' object has no attribute 'encode''

AttributeError                            Traceback (most recent call last)
<ipython-input-16-80c4e8c0de28> in <module>()
      1 sql = 'Select * FROM table'
----> 2 GeoDF= gpd.GeoDataFrame.from_postgis(sql, engine)

~\AppData\Local\Continuum\anaconda3\lib\site-packages\geopandas\geodataframe.py in from_postgis(cls, sql, con, geom_col, crs, index_col, coerce_float, params)
    235         """
    236         return geopandas.io.sql.read_postgis(sql, con, geom_col, crs, index_col,
--> 237                      coerce_float, params)
    238
    239     def to_json(self, na='null', show_bbox=False, **kwargs):

~\AppData\Local\Continuum\anaconda3\lib\site-packages\geopandas\io\sql.py in read_postgis(sql, con, geom_col, crs, index_col, coerce_float, params)
     40     wkb_geoms = df[geom_col]
     41
---> 42     s = wkb_geoms.apply(lambda x: shapely.wkb.loads(binascii.unhexlify(x.encode())))
     43
     44     df[geom_col] = GeoSeries(s)

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pandas\core\series.py in apply(self, func, convert_dtype, args, **kwds)
   3190             else:
   3191                 values = self.astype(object).values
-> 3192                 mapped = lib.map_infer(values, f, convert=convert_dtype)
   3193
   3194         if len(mapped) and isinstance(mapped[0], Series):

pandas/_libs/src\inference.pyx in pandas._libs.lib.map_infer()

~\AppData\Local\Continuum\anaconda3\lib\site-packages\geopandas\io\sql.py in <lambda>(x)
     40     wkb_geoms = df[geom_col]
     41
---> 42     s = wkb_geoms.apply(lambda x: shapely.wkb.loads(binascii.unhexlify(x.encode())))
     43
     44     df[geom_col] = GeoSeries(s)

AttributeError: 'NoneType' object has no attribute 'encode'


texts: QGIS virtual raster - insufficient permissions

texts:
I'm trying to create a virtual raster in QGIS 2.6.1 64b on Windows, but I get the message:

  "The process failed to start. Either the invoked program is missing,
  or you may have insufficient permissions to invoke the program."

I try to read in a whole directory, selecting also to "Recurse subdirectories"
it works for a part of my datasets, but not for all.
I thought it is related to the the size of the datasets, so I tried dividing the bigger datasets in smaller directories. but it did not work.


texts: How to buffer a line with different distances for left/right sides

texts: How can I buffer a line in QGIS using one attribute for the left side of the line and one attribute for the right side of the line? I am trying to buffer a river with its bank width, which is different on the right and left side.


texts:
Good day all,

I've successfully used the Qgis2Threejs plugin to view a 3D scene (pretty basic but it does the job).

Does anyone know if it's possible to create a fly-through of the model? Something like this; https://youtu.be/PbDWSKPBXNI


texts: Monitoring Progress of projectRaster in R

texts:
I am new to R and I am having a few issues with this piece of code. The goal of it is to take set of daily temperature rasters of the United States, reproject them, and clip them. At the moment, I am doing this by creating a raster stack, projecting it, and then clipping it. When I run this block of code, I don't get any errors, but it gets hung up on the projectRaster step (at the moment I am trying to reproject a stack of 365 days).

library(raster)
library(rgdal)
#import area of interest file
boundary <- readOGR(dsn = 'D:/Scratch/DataConvert/Boundary' , layer =
'AreaOfInterest' )

#define spatial reference
sr <- ('+proj=lcc +lat_1=45 +lat_2=49 +lat_0=44.25 +lon_0=-109.5
+x_0=600000.0000000001 +y_0=0 +ellps=GRS80 +datum=NAD83
+to_meter=0.3048006096012192 +no_defs ')

#project boundary
pr <- spTransform(boundary, crs(sr))

#Create 'images' list that is composed of daily temps for the entire US
images <- list.files(path = 'D:/Scratch/DataConvert/Daily/', pattern =
'*.tif' , full.names = TRUE, recursive = FALSE)

#Create a raster stack of the images list
s <- stack(images)

#Project the raster stack
s.proj <- projectRaster(s, crs = sr)

#Crop the raster stack to the area of interest
s.crop <- crop(s.proj, extent(pr))

#Save a clipped, projected image for each day
writeRaster(s.crop, filename = 'D:/Scratch/DataConvert/JAF_Clipped/clip',
        bylayer =TRUE, progress='text' , suffix=names(s), format='GTiff' ,
        overwrite=TRUE)

It takes an extremely long time (>24hours) for this segment to run and I want to be able to monitor its progress. I read through this article that describes setting up a progress bar but I can't seem to make it work. This was one of the many ways I tried to implement the progress bar.

total <- 100
pb <- txtProgressBar(min = 0, max = total, style = 3)
s.proj <- projectRaster(s, crs = sr)
for(i in 1:total){
   Sys.sleep(0.1)
   # update progress bar
   setTxtProgressBar(pb, i)
}
close(pb)

Does anyone have any insight into incorporating a progress bar with the projectRaster function? Or does anyone know of a more efficient way to complete this task?


texts:
I'm trying to snap a line theme to a point theme using PostGIS with a filter so they only snap to records with the same reference number.  The query I'm using is as follows (this is a sample of the data and without the query, just trying to get the snapping working!)

CREATE TABLE public.snappedfeatures3mV13_TABLE AS
SELECT
    st_snap(featuresshape_cut.geom, "nodes_cut ".geom, 2::double  precision) AS geom,
    featuresshape_cut.id ,
    featuresshape_cut.site_code
FROM featuresshape_cut
JOIN "nodes_cut " ON featuresshape_cut.site_code::text = "nodes_cut ".site_code::text

No errors received, seems to generate the output file fine, however the output is multiple features linking a snapped point and the two closest points, image below.  Green dots are the features snapped too, the black line the original and red is the output with arrows to indicate ends of line segments.  Its almost like its calculating the snaps it needs to do then stopping?

I have

Converted the original features to be snapped to from lines to points - no effect
Cut the dataset down to a minimum
Checked all inputs have valid geometry
Generated new primary keys for all inputs
played with snapping tolerance to ensure snapping is working as expected
Checked all CRS are the same

No errors in PostGIS.  I'm obviously missing something, what is it?

Demonstration of a single feature from the new file, there are 3 nodes in each record, the "middle" has snapped to the line as expected and the start and end are on nodes in the original line.


texts:
I created overlapping geometries and I would like to know the area of the intersection.

import geopandas as gpd
from shapely.geometry import Polygon

polys1 = gpd.GeoSeries([Polygon([(0,0), (0,1), (1,1), (1,0)]),
                        Polygon([(89,89), (89,90), (90,90), (90,89)]),

                        ])

polys2 = gpd.GeoSeries([Polygon([(0.5,0.5), (0.5,1.5), (1.5,1.5), (1.5,0.5)]),
                        Polygon([(88.5,88.5), (88.5,89.5), (89.5,89.5), (89.5,88.5)]),
                        ])

df1 = gpd.GeoDataFrame({'geometry': polys1, 'seed_value':[10,10]})
df2 = gpd.GeoDataFrame({'geometry': polys2})

merged_df = gpd.overlay(df2, df1, how='intersection')

print(merged_df['geometry'][0].area)
print(merged_df['geometry'][1].area)

One intersection is located at the equator (left) and the other at a pole (right). The code returns an area for both intersections of 0.25.

Shouldn't the area be different because of the different length of arcminutes at the equator and at the poles? How should I interpret the result and how can I compute the 'real' area?

If I would like to pass on values from one geometry to another geometry weighted by the overlapping area, what values do I have to compute to take into account the different lengths of arc-minutes at the poles and at the equator?


texts: CartoCSS: labels along inside of polygon boundary

texts:
I'm using CartoCSS to render maps that include large polygons (like national parks). At higher zoom levels I want labels along the inside of the polygon boundaries, like in OpenStreetMap:

I reviewed the openstreetmap-carto stylesheets but could not find how this is done. Using text-placement: line just places text on the line (this is how streets are labelled). I need it inside the line. Any ideas?


texts: Filter out the least cloudy images in Sentinel Google Earth Engine

texts:
I am trying to filter out the least cloudy Sentinal 2 images and display them on as per the AOI boundary. I am working on the following code, but it is not filtering out the least cloudy images or any method to remove the cloud from images. I had gone through several posts but not found any solution yet.

var shp=ee.FeatureCollection('users/Desktop/AOI');
Map.addLayer(shp);
var sentinel=ee.ImageCollection('COPERNICUS/S2_SR').filterDate('2019-12-10','2019-12-30')
.select(['B8','B4','B3','B2']).filterBounds(shp).sort('CLOUDY_PIXEL_PERCENTAGE', false);
var s2a=sentinel.mosaic().clip(shp);
print(sentinel);
//var clipped = sentinel.map(cloudmask).map(clipper);
//var leastCloudy = ee.Image(sentinel.sort('CLOUDY_PIXEL_PERCENTAGE').first());
Map.addLayer(s2a,{min:0, max:4000}, 'Output');


texts:
I'm working with some raster data I downloaded. The problem is I can't manage to find the correct projection for this file. It should work with WGS84 (as shown in the image below) but it's displaced towards the south relative to the borderline. The QGIS plugin "Projestions" doesn't suggest any either.
Any ideas about how this can be solved?

The output of gdalifno:
Warning 1: Latitude grid not spaced evenly.  Setting projection for grid spacing is within 0.1 degrees threshold.Driver: netCDF/Network Common Data Format
Size is 631, 539
Origin = (-5.811904761904762,44.722040824551812)
Pixel Size = (0.023809523809524,-0.018239630129158)
Metadata:
  lat_rho#axis=Y
  lat_rho#field=lat_rho, scalar
  lat_rho#long_name=latitude of RHO-points
  lat_rho#standard_name=latitude
  lat_rho#units=degree_north
  lat_uv#axis=Y
  lat_uv#field=lat_uv, scalar
  lat_uv#long_name=latitude of UV-points
  lat_uv#standard_name=latitude
  lat_uv#units=degree_north
  lon_rho#axis=X
  lon_rho#field=lon_rho, scalar
  lon_rho#long_name=longitude of RHO-points
  lon_rho#standard_name=longitude
  lon_rho#units=degree_east
  lon_uv#axis=X
  lon_uv#field=lon_uv, scalar
  lon_uv#long_name=longitude of UV-points
  lon_uv#standard_name=longitude
  lon_uv#units=degree_east
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#date_created=07-Aug-2020 10:06:17 CET
  NC_GLOBAL#file=WMOP surface his file
  NC_GLOBAL#format=netCDF-3 classic file
  NC_GLOBAL#fullname=Balearic Islands Coastal Observing and Forecasting System
  NC_GLOBAL#institution=SOCIB
  NC_GLOBAL#publisher_email=info@socib.es
  NC_GLOBAL#publisher_name=SOCIB
  NC_GLOBAL#publisher_url=www.socib.es
  NC_GLOBAL#version=WMOPv2.2r ROMS Forecast (with data assimilation)
  NETCDF_DIM_EXTRA={ocean_time}
  NETCDF_DIM_ocean_time_DEF={25,6}
  NETCDF_DIM_ocean_time_VALUES={1647561600,1647572400,1647583200,1647594000,1647604800,1647615600,1647626400,1647637200,1647648000,1647658800,1647669600,1647680400,1647691200,1647702000,1647712800,1647723600,1647734400,1647745200,1647756000,1647766800,1647777600,1647788400,1647799200,1647810000,1647820800}
  ocean_time#axis=T
  ocean_time#calendar=gregorian
  ocean_time#field=time, scalar, series
  ocean_time#long_name=time since initialization
  ocean_time#units=seconds since 1968-05-23 00:00:00 GMT
  v#coordinates=ocean_time lat_uv lon_uv
  v#field=v-velocity, scalar, series
  v#long_name=v-momentum component
  v#standard_name=northward_sea_surface_velocity
  v#time=ocean_time
  v#units=meter second-1
  v#_FillValue=-nan
Geolocation:
  LINE_OFFSET=0
  LINE_STEP=1
  PIXEL_OFFSET=0
  PIXEL_STEP=1
  SRS=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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  X_BAND=1
  X_DATASET=NETCDF:/roms_wmop_surface_20200807.nc":lon_uv
  Y_BAND=1
  Y_DATASET=NETCDF:/roms_wmop_surface_20200807.nc":lat_uv
Corner Coordinates:
Upper Left  (  -5.8119048,  44.7220408)
Lower Left  (  -5.8119048,  34.8908802)
Upper Right (   9.2119048,  44.7220408)
Lower Right (   9.2119048,  34.8908802)
Center      (   1.7000000,  39.8064605)
Band 1 Block=631x1 Type=Float32, ColorInterp=Undefined
  Min=-0.797 Max=0.995
  Minimum=-0.797, Maximum=0.995, Mean=0.036, StdDev=0.180
  NoData Value=nan
  Unit Type: meter second-1
  Metadata:
    coordinates=ocean_time lat_uv lon_uv
    field=v-velocity, scalar, series
    long_name=v-momentum component
    NETCDF_DIM_ocean_time=1647561600
    NETCDF_VARNAME=v
    standard_name=northward_sea_surface_velocity
    STATISTICS_MAXIMUM=0.99483668804169
    STATISTICS_MEAN=0.036012467945364
    STATISTICS_MINIMUM=-0.79715889692307
    STATISTICS_STDDEV=0.18019412315803
    time=ocean_time
    units=meter second-1
    _FillValue=-nan


texts: Calculating slope of numpy array using gdal.DEMProcessing

texts:
I am trying to compute the slope of a numpy array using gdal.DEMProcessing. Here is an example of what I am trying. Here in the example I am loading a portion of a DEM file, when I try to compute the slope I get out an array of -9999.

dem_ds = gdal.Open(str(dem_filename))
dem_arr = dem_ds.ReadAsArray(xoff=3584, yoff=0, xsize=512, ysize=512)

dst_filename = '/vsimem/slope.tif'

tmp_ds = gdal_array.OpenArray(dem_arr)
_ = gdal.DEMProcessing(dst_filename, tmp_ds, 'slope')
slope = gdal.Open(dst_filename).ReadAsArray()
slope.mean()

  -9999.0

It works by passing dem_ds instead of tmp_ds however the point would be to compute the slope of a general numpy array (without georeference) and not the full dataset as loaded from file.


texts:
I am trying to enter in my maps application (GeoServer, PostGIS, OpenLayers) but it is empty and I get an error with Firebug: TypeError: OpenLayers.Map is not a constructor in my map.js script in line:

var map = new OpenLayers.Map({
    projection: "EPSG:900913",
    displayProjection: "EPSG:4326",
    numZoomLevels: 7});`

What is wrong?

EDIT:
This is my header html

<script type="text/javascript" src="http://extjs.cachefly.net/ext-3.4.0/adapter/ext/ext-base.js"></script>
<script type="text/javascript" src="http://extjs.cachefly.net/ext-3.4.0/ext-all.js"></script>
<link rel="stylesheet" type="text/css" href="http://extjs.cachefly.net/ext-3.4.0/resources/css/ext-all.css" />
<link rel="stylesheet" type="text/css" href="http://extjs.cachefly.net/ext-3.4.0/examples/shared/examples.css" />
<script src="http://openlayers.org/dev/OpenLayers.js"></script>
<script type="text/javascript" src="./script/GeoExt.js"></script>
<script src="http://maps.google.com/maps/api/js?v=3&amp;sensor=false"></script>

<link rel="stylesheet" href="http://openlayers.org/dev/theme/default/style.css" type="text/css">
<link rel="stylesheet" href="http://openlayers.org/dev/theme/default/google.css" type="text/css">

<script type="text/javascript" src="map.js"></script>

The problem occurred suddenly while I looking the layers before added. I clicked in option Reload feature type in a particular layer page, Could this have something to do?


texts: How to detect if a point lies inside a building?

texts:
My main motive is to detect if there's in any infrastructure nearby to a given point in a certain radius. For this I am using Places api and using nearbySearch method of it. The problem in this is that it only detects any infrastructure that is already marked by someone on Google Map's.
Eg: http://web.iiit.ac.in/~jayesh.lahori/map.png , blue point is itself inside the building, this will not be detected because nothing has been marked at this spot.

Any suggestions ??


texts: Auto Adding Non Spatial Table into QGIS template from postgres PGAdmin (customize python plugin)

texts:
I am working with customize plugin python Qgis (version 2.18). This plugin login from QGIS dashboard and calling spatial data from postgresql PGAdmin (version 4 v4).
The moment user key in username and password in QGIS customize plugin, all spatial layer will be auto sort in layer panel accordingly.
My prob comes on how to auto call non-spatial layer into the layer panel. I make it a same query with spatial layer in postgres;
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('horizontal_weephole', 'Horizontal Weephole', 'CERUN', '1', 't');
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('instrument_cerun', 'Instrument Cerun', 'CERUN', '2', 't');
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('sistem_pengukuhan', 'Sistem Pengukuhan', 'CERUN', '3', 't');
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('sistem_perparitan_cerun', 'Sistem Perparitan Cerun', 'CERUN', '4', 't');
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('slope', 'Slope', 'CERUN', '5', 't');

INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('slopeinfo', 'Slope Info', 'CERUN', '6', 't');
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('services_facilities', 'Services Facilities', 'CERUN', '7', 't');
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('sign_of_distress', 'Sign of Distress', 'CERUN', '8', 't');
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('surface_subsurface', 'Surface Subsurface', 'CERUN', '9', 't');
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('engineering_judgement', 'Engineering Judgement', 'CERUN', '10', 't');
INSERT INTO "public"."mc_kej_cer_v_layers" (layer_name, layer_alias, parent_group, order_in_parent, visibility_at_start) VALUES ('structure', 'Structure', 'CERUN', '11', 't');

and also grant the same permission as spatial layer. But the result comes out with no table appear in the layer panel.

Supposedly it will be appear like this;

The second picture get from manually added the tables.
After that, all the tables need to join with spatial layer Slope using 'slopeid' as a primary key. All this can be done manually and save the style. But the problem is it become temporary work. When user log out from the customize application and login back, all join setting and the non spatial table disappear.
Is it possible if join table with spatial table will be done in PostgreSQL first and produce the new spatial table?
Then calling the new join table slope into QGIS layer panel.
The new table can be 'View Update Delete  attribute?
Can you share me the query code?


texts: 3D wkt geometries - Parse in .Net suggests alternate form

texts:
I am having some issues with a built in .net framework class (System.Data.Spatial.DbGeometry in assembly System.Data.Entity). It turns out that this class cannot parse 3D geometries when they are denoted using the correct (at least to my understanding) form.

According to Section 7.2.3 of OpenGIS Implementation Specification for Geographic information - Simple feature access - Part 1: Common architecture

The correct representation of a 3D point is: "Point Z (10 10 5)". However .Net doesn't understand the 'Z'. If I omit the 'Z' and try to parse "Point (10 10 5)", I get a 3d point.

My question is, is there such an alternative WKT-specification that allows omitting the 'Z' (same goes for measures 'M')?


texts: Is there a way of doing ST_SnapToGrid in MySql?

texts:
Given that MySQL 5.7 has spatial data types
Is there an equivalent to PostGIS's ST_SnapToGrid in MySQL?

Looking at the reference:

  http://dev.mysql.com/doc/refman/5.7/en/spatial-function-reference.html

I can't find one. Is there something else I could use that would do the same thing?


texts: I have 0 experience in AutoCAD and need to convert some layers to .shp I can use and store in ArcMap. Would anyone be willing to give a step-by-step procedure on this process?


texts: Migrating Metadata Storage from SQL Server to ArcGIS Online

texts:
Currently we have somewhere around 2500 metadata records in a SQL Server table to relate to an online mapping tool. I am trying to formulate some sort of idea of how I can take all those records out and import them into the metadata standard using ArcGIS Online.

Ideally, this would be solved by extracting the SQL Server table records as an XML for each record, and then rename or map the SQL Server fields to the applicable AGS Online metadata fields. Lastly, the new XML would be imported somehow into AGS Online, populating the relevant fields defined before.

Since there's a lot of records, I would like to automate at least part of this. Does anyone have a medium or high level idea of how this could happen?


texts: Making grid files (minimum curvature) from point data in QGIS?

texts:
Is there any way to make a grid file from point data (shape file) in QGIS?

I mean by grid a colored raster file like Geosoft .grd or ERMAPPER .ers file or like.

Minimum curvature would be especially nice to have a gridding algorithm.


texts: How to create a Qgis python script to manipulate fields with some parameters?

texts:
I'm just beginning with QGIS.

I have a shape file with 10 fields field1, ..., field10, and many thousands of data.

I will like to automatize the next sequence (in a python script or a model):

load the shape file
create a field field_add defined as field1 + field2
value in field3 is Text1, Text2, Text3 or Text4 (if Text1, coeff=a; if Text2, coeff=b; if Text3, coeff=c; if Text4, coeff=d; else coeff=e)
create a field field_age defined as field4 * coeff
create a field field_result defined as field_age^m * field_add^n
create a field field_class defined as 1 if field_result < floor1, 3 if field_result > floor2, else 2.

After that, I will visualize each polygon with its field_class color (white for 1, green for 2 and red for 3).

All the coefficients (a,b, c, d, e, m, n, floor1 and floor2) have to be changed quickly, to update the chart.

I insert part of the table:

Please, could anyone indicate to me the best way to proceed?


texts: Is there any public ArcGIS Server that uses HTTP authentication that can be used for testing?

texts:
I have the requirement to connect to an ArcGIS service that use Http authentication, currently I do not have an EDN or Server license to setup an environment like that.
Are there any public server that can provide this services with some kind of free user?


texts: Creating a raster surface of overlaps from vector polygon shapefiles?

texts:
Basically what I'm trying to do (in ArcGIS 10) is to compare three vector files in order to create a 'comparison' raster layer. I tried using Data Comparison (under Data Management tools) but it didn't output what I was looking for. I basically wanted to create a raster surface that shows where each type of polygon overlaps. So for example if  layer a and b overlap at some point, the raster pixel value will reflect that (and if a, b and c overlap, etc.). I can't figure it out and I've been browsing ArcMap's help for ages.

Thank you!


texts: Producing raster mosaic dataset using spatial query on another mosaic in ArcGIS for Desktop?

texts:
Using ArcGIS 10.1 (Build 3143) I would like to create a relatively small raster mosaic dataset (the size of a city) based on a larger, existing one (the size of a state) using a single feature class polygon as a spatial query.  Is this possible?

Clipping the raster mosaic with Clip (Data Management Tools > Raster > Raster Processing > Clip) or Extract by Mask (Spatial Analyst Tools > Extraction > Extract by Mask) produce low-resolution rasters instead of mosaics.

Over the past few weeks I think I remember seeing the option to perform the operation described above, but now I cannot find it.  Ultimately, I will publish the raster mosaic dataset in ArcGIS Server.


texts: Changing default zoom in Google embedded map?

texts:
I am having difficulties to change the default zoom level when I open the iframe below in my site. Can anyone please help me to set zoom level when following iframe is opened?

<iframe src="http://mapsengine.google.com/map/u/0/embed?mid=zYk_BmXXTdxw.k83snDbX5SJc" width="745" height="480"></iframe>

This is the newest version of Google Maps.


texts: Is it possibleto bar graph rainfall data in parallel with with a temporal animation

texts: I have an animation of red squirrel sightings in Clocaenog Forest over time (one year).  I would like to display a bar graph of rainfall and temperature for each day as the animation displays.  I know there are graphing facilities in QGIS, can they be displayed in a temporal animation?


texts: Finding the Python script used by a custom toolbar button in ArcMap

texts:
I've been asked by somebody at work to design a toolbar add-in that works just the same as an existing one, but with a different set of data. However, the person who wrote the original add-in for him has left the company, and did not document his work at all.

The only thing available to me is the MXD. How do I find what Python script(s) is triggered by a toolbar add-in, and where the script lives? (The requester is working in ArcMap 9.3, and while I have ArcMap 10.2 available to me, I have to design to 9.3 so he'll be able to use it.)

An obvious workaround is just writing something new from scratch. The desired workflow is fairly simple: iterate through a bunch of groups of data and make a PDF, thereby automating the process of making 100+ maps. But just like this analyst is simplifying his job by using a one-click button, I'd like to simplify my job by basing the new script on the existing, proven script :-)


texts: ESRI JS API hosted locally but requests still go to https://js.arcgis.com?

texts:
In my current project we have hosted the ESRI JS API within the organisation's IIS.

For most of the part the deployed application works fine, but we notice that certain calls are being made to https://js.arcgis.com despite the corresponding JS files being present in the local deployment.

I have deployed the ESRI API as per instructions within their readme/install documentation. I have not completed the steps under "Change the ArcGIS Services Directory "View In JavaScript" URL" section, as I think these are optional and only apply when viewing the deployed map services from within the REST services directory.

Any ideas why the external calls are being made?


texts: calculate area of all polygons in shapefile in arcgis

texts: I have shapefile containing lot of polygons, I want to know area of all of them together. Do you please know about some few clicks tool for my problem?


texts:
I'm quite new to PostGIS and I'm trying to understand how I can optimise a simple geo query that seems to be a bit slow.

I'm using GeoDjango, but I've tried to perform the queries directly on PostgreSQL using pgAdmin, with the same results in terms of performance.

I've a table called busstop which main geo field is a Point and I've a location table which main field is a MultiPolygon.

Given a point (which I know is within 10-15 Km range from the other records) I'm trying these queries in GeoDjango (query: "give me all the records within 2Km from the base point")

Location.objects.filter(geom__distance_lt=(base_point, D(Km=2))).count()

and

BusStop.objects.filter(point__distance_lt=(base_point, D(Km=2))).count()

Just to give you a reference, busstop contains ~18000 records while location contains ~14000. The average time to compute these queries was ~4 seconds.

What I've tried to do:

I've recreated the location geom field using this:

geom = models.MultiPolygonField(geography=True, spatial_index=True)
I've used pgAdmin -> Manage to vacuum, analyze, reindex both location and busstop table.

Result: query time have been reduced from ~4 seconds to ~2 seconds. It's something good, but it's still not the performance I would expect for so little data.

I'm pretty sure I'm doing something wrong or I'm missing something trivial that would let me have better performance. Has anyone any idea about what else I could do to improve these simple queries time?

I've tried to drop the database, re-create it, re-import the initial data and executed the queries again, with this results:

busstop query time went from 2 seconds to 183 ms
location query time went from 2 seconds to 1 second

I still think the location query could be further optimised.

This is the query executed for location:

SELECT "api_location"."id", "api_location"."name", "api_location"."point", "api_location"."geom", "api_location"."authority", "api_location"."owner", "api_location"."uprn", "api_location"."unique_asset_id" FROM "api_location" WHERE ST_Distance("api_location"."geom", ST_GeomFromEWKB('\001\001\000\000 \346\020\000\000\210\215\325\247\202\367\001\300q\216\350\015\272\301J@'::bytea)) < 2000.0

Result of explain is:

Seq Scan on api_location  (cost=0.00..4766.05 rows=3815 width=1942)
Filter: (_st_distance(geom, '0101000020E6100000888DD5A782F701C0718EE80DBAC14A40'::geography, '0'::double precision, true) < '2000'::double precision)


texts:
I'm trying to rasterize a polygon layer which has two columns : id and code.
I would like each polygon to get value according to the code field.
I have used the tool rasterize from the raster menu

but it keeps failing with the next error message:

Loading resulting layers The following layers were not correctly
generated.C:/Users/hello/this/is/Tmy/processing/fmoreNumbers/OUTPUT.tifYou
can check the 'Log Messages Panel' in QGIS main window to find more
information about the execution of the algorithm.

why is it failing all the time? what can be the problem? and can I rasterize also based on categorical values?
Edit: log message:
Algorithm 'Rasterize (vector to raster)' starting…
Input parameters:
{ 'BURN' : 0, 'DATA_TYPE' : 5, 'EXTENT' : '810660.0,844250.0,8123830.0,8145680.0 [EPSG:32721]', 'FIELD' : 'class', 'HEIGHT' : 0, 'INIT' : None, 'INPUT' : 'MultiPolygon?crs=EPSG:4326&field=Plot:string(254,0)&field=Crop:string(254,0)&field=Start:string(254,0)&field=End:string(254,0)&field=Grower:string(254,0)&field=Group:string(254,0)&field=Region:string(254,0)&field=Area%2C%20ha:double(23,15)&field=Season%20yea:string(254,0)&field=Code:string(254,0)&field=class:long(10,0)&uid={8ab7aaf2-0e7c-47f2-82d6-dd2abdb6d30c}', 'INVERT' : False, 'NODATA' : 0, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'UNITS' : 0, 'WIDTH' : 0 }
GDAL command:
gdal_rasterize -l INPUT -a class -ts 0.0 0.0 -a_nodata 0.0 -te 810660.0 8123830.0 844250.0 8145680.0 -ot Float32 -of GTiff C:/my/path/INPUT.shp C:]my/path/OUTPUT.tif
GDAL command output:
ERROR 1: Wrong value for -outsize parameter.
Usage: gdal_rasterize [-b band]* [-i] [-at]
{[-burn value]* | [-a attribute_name] | [-3d]} [-add]
[-l layername]* [-where expression] [-sql select_statement]
[-dialect dialect] [-of format] [-a_srs srs_def] [-to "NAME=VALUE"]*
[-co "NAME=VALUE"]* [-a_nodata value] [-init value]*
[-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]
[-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
CInt16/CInt32/CFloat32/CFloat64}] [-optim {[AUTO]/VECTOR/RASTER}] [-q]
<src_datasource> <dst_filename>
Execution completed in 0.33 seconds


texts: Which VS projects are needed for a custom QGIS build?

texts:
Environment:

VS 2010
Win 7
QT 4.8.6
QGIS git Src
OSgeoW4

I have built the source for QGIS, in order to use it as a base for a custom application.
I wanted to know what projects are 100% needed, after the first build, just to make changes to the ui?

212 projects loaded in VS is a lot of memory, so any guidance would be great.


texts: Getting house neighbors using GeoDjango?

texts:
i have a house model with these properties

centerpoint = PointField(srid=0, null=True, blank=True)
front_left = PointField(srid=0, null=True, blank=True)
front_right = PointField(srid=0, null=True, blank=True)
rear_left = PointField(srid=0, null=True, blank=True)
rear_right = PointField(srid=0, null=True, blank=True)
polygon = PolygonField(srid=0, null=True, blank=True)
direction = FloatField(null=True, blank=True)  # the direction this house faces, in degrees

and i want to get the houses in front of a specific house... im having a hard time trying to use Distance and polygon__dwithin

i need to get the ones that are touching and the ones that have relative direction of less than 160 degrees:

for that i do:

House.objects.filter(polygon__touches=house_instance.polygon, direction__lt=160)

how can i modify this query to get the houses that are in front maximum one block and the ones that are one house away from my reference house?

also is it possible to get the neighbors by using my polygon field to get the closest polygon fields from it?


texts: Identifying bands of satellite images qgis

texts:
I have a series of satellite images and am not sure which satellite took them.
I know there are 8 bands but not which band relates to which spectral color. Here is the metadata of one image :

Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
/home/kai/sat images from digital globe/054892521020_01 katchall 15,13 p,m/054892521020_01_P002_MUL/13JAN07045052-M2AS_R2C1-054892521020_01_P002.TIF
AREA_OR_POINT=Area
TIFFTAG_COPYRIGHT=(C) COPYRIGHT 2015 DigitalGlobe, Inc., Longmont CO USA 80503
TIFFTAG_DATETIME=2015:12:11 16:13:54
TIFFTAG_IMAGEDESCRIPTION={ bandList = [ 6; 2; 3; 7; 4; 8; 5; 9; ] }
TIFFTAG_MAXSAMPLEVALUE=2047
TIFFTAG_MINSAMPLEVALUE=0
Band 1
STATISTICS_MAXIMUM=1478
STATISTICS_MEAN=257,95138713292
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=157,66542923619
Band 2
STATISTICS_MAXIMUM=1764
STATISTICS_MEAN=202,60800252642
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=128,2493741809
Band 3
STATISTICS_MAXIMUM=2047
STATISTICS_MEAN=201,06180054801
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=146,04955395515
Band 4
STATISTICS_MAXIMUM=2047
STATISTICS_MEAN=126,90268325806
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=112,31556464568
Band 5
Band 6
Band 7
Band 8
Dimensions
X: 4096 Y: 4096 Bands: 8
Origin
531286,877766
Pixel Size
2,-2
No Data Value
*NoDataValue not set*
Data Type
UInt16 - Sixteen bit unsigned integer
Pyramid overviews
Layer Spatial Reference System
+proj=utm +zone=46 +datum=WGS84 +units=m +no_defs
Layer Extent (layer original source projection)
531286.0000000000000000,869574.0000000000000000 : 539478.0000000000000000,877766.0000000000000000
Band
Band 1
Band No
1
No Stats
No stats collected yet
Band
Band 2
Band No
2
No Stats
No stats collected yet
Band
Band 3
Band No
3
No Stats
No stats collected yet
Band
Band 4
Band No
4
No Stats
No stats collected yet
Band
Band 5
Band No
5
No Stats
No stats collected yet
Band
Band 6
Band No
6
No Stats
No stats collected yet
Band
Band 7
Band No
7
No Stats
No stats collected yet
Band
Band 8
Band No
8
No Stats
No stats collected yet


texts:
I was given some data sets to create maps in QGIS. Most of them were shape files, but some of them are in a format I am not sure about. For example it is a folder called "Rivers", containing AAT.dbf, abs, abs, arc, ARF, ARX, BND.dbf, CNT, CNX, LAB, MSK, NRF and other files. Are these ArcInfo files? I am attaching a screenshot.

Could somebody help me how to open them in QGIS?


texts: Adding Parameter Descriptions to Script Tool

texts:
I've made a simple tool that accepts three parameters: an input file (feature class), an output folder (folder), and an output name (string). These are all required input with no filters. The tool's main description correctly displays on the right pane of the tool window in Arc, but the input parameter descriptions are not showing up when their fields are clicked on. I tried editing the Item Description in ArcCatalog, but this only brings up an empty window with nothing to edit. What am I missing here?

This is the tool's

The tool executes successfully.

Ideally I'd like to be able to access/edit documentation completely within the script, but I don't know how to start there. This is all saved within a .xml file somewhere, right? If so, where is that stored?

EDIT: I'm using ArcGIS 10.2.2 for Desktop. Also, in case it's relevant, I made this tool by making a new toolbox (.tbx), right-click -> Add -> Script... And I browsed to my .py file.

#!usr\bin\python

import arcpy, os, csv

# Read the parameter values:
infile = arcpy.GetParameterAsText(0)
outfolder = arcpy.GetParameterAsText(1)
outname = arcpy.GetParameterAsText(2)

if outname[:-4] != '.csv':
    outname += '.csv'
# etc...


texts: How to design a database for continents, countries, regions, cities and POIs?

texts:
I'm brand new to GIS programming and I am designing a GIS application. Target is to create system with continents, countries, regions (including states, sub-regions, provinces), cities and places in cities. Each of this elements will contain some text information and related stuff. As database we are going to use PostgreSQL with PostGIS.

My question is how to design database for this system? I was thinking of 2 tables polygons and points, but I'm not sure if it's good way of thinking.

What we need is some hierarchical base for relationship between countries, cities, regions, etc.

Main goal of the application will be collecting tourist data from many sources about specific cities, regions, countries and so on, and store it in database. Let's say we will have city Rome, we collect some info about this city into database from couple sources. And we need to know that Rome is in province Rome, sub-region Lazio in region Lazio, country Italy. So system should be flexible to allow create such relation from real world.