texts:
I have about 2 massive layer files with polygons and like 10 layer files with points that I would like to "Merge" together separately (polygons in one and points in another) but for some reason It removes my firs column with important ID's and replace it with randomly generated ObjID? How can I prevent this from happening? Another problem that I found was that everytime I tried to marge the 2 massive layers with points it finished the merging prosess but when I look at the data in the attribute table it always stopped at like 6025 or something. I think it should be more like ca 10000.

I want to do this as I currently have multiple layer files that I want to run through both "generate near table" and "make feature layer". The issue is that "make feature layer" can only do 1 layer. Hence I thought It might work if I can just have all the polygons layers as 1 superlayer and all the points as another super layer?


texts: ArcMap Using Python to Rename a Feature Class based on Print Results

texts:
I would like to rename a feature class based on the EDITOR and DATE fields. I wrote this code:

    import arcpy

fc = "C:\Users\PythonTesting.gdb\MyFile"
f1, f2, f3 = "EDITOR", "DATE", "OBJECTID"
clause = arcpy.AddFieldDelimiters(fc, f3) + "= 1"
for row in sorted (arcpy.da.SearchCursor (fc, [f1, f2], clause)):
     print ("{0}_{1}".format(row[0], row[1]))

And it works great, prints out the one line of results I want. But now I can't figure out how to use those results as the output for a rename. I have used the rename before but with out such an expression.

I checked out this one thinking it would help but I still can't figure it out.
How to rename feature classes based on field values within the attribute table?


texts:
I have some data that contains crab males population based on latitude and longitude. I am trying to plot the population distribution on map (not point density). For example I like to show which areas are more populated and where less.
Here is my data

     lon    lat NUM_MALES_KEPT
1 144.35 -41.51              4
2 144.35 -41.52             22
3 144.39 -41.56              3
4 144.38 -41.53              5
5 144.38 -41.52              3
6 144.37 -41.54             10

I used geom_density code but this code maps data based on point density not population density.


texts: Running script over every file in folder using GDAL in Python?

texts:
I'm using GDAL to calculate NDVI mean and max and to find other information about .tif Landsat images acquired from Earth Explorer.  I stacked the bands into a single .tif image and have a working python script to find everything I want to know about the image, however I have many .tif files and I don't want to run the script individually for each file, and I want to make my script automatically run the code over every .tif file I have in the same folder.

The beginning of my script, and how I am accessing the single .tif file looks like this:

from __future__ import print_function

from osgeo import gdal

print("GDAL's version is: " + gdal.__version__)

print (gdal)

print(gdal.GDT_Byte)

dataset = gdal.Open('/Path/To/Folder/Containing/.tif', gdal.GA_ReadOnly)

How would I get Python to run the script over every file in this folder automatically?


texts: Including elevation values in polyline (to export it as ASCII table)?

texts:
I'm using a code to simulate flood events and I need to build two main files: one with the river channel width and another one with the bed elevation (both in ARC-info ascii raster format).

I've managed to create a new field in the attributes table and insert a constant width value, but when it comes to generate the raster grid with the bed elevation I'm lost. I've tried "Add surface information", "join" and "interpolate" tools so far, but it doesn't not work as I expected.

The final outcome has to be a ascii file with NoData cells (=-9999) for those areas where the river channel does not exist and bed elevation values for those cells where the river is present.

I have beginner-intermediate programming knowledge (VBA, DOS and a bit of Python mainly) and I'm working with ArcGIS 10.4.1; the final grid has to contain 1500 columns and 1000 rows, since I'm working with 6 tiles of 2 m resolution DEM data.


texts: point feature to raster with count of unique points

texts:
I am using ArcGIS 10.1 with Python 2.7
I have a point layer with thousands of points that each hold an id value of its owner. One owner can have multiple points, meaning that points can have the same owner_id.
I want to create a raster that counts the number of unique owner_ids in each cell.

Is there a way to create a raster with a value containing the count of unique owner_ids?
In other words, I want a raster with the number of unique owners instead of a raster with the number of points in each cell.

I have tried creating separate raster grids for each individual owner and its points (found it somewhere on another forum), but since my data in total is more than a million different owner_ids, this will take way too long.


texts:
I got credential by using this command.

earthengine authenticate

I have linked with my account and saved authentication token successfully. I can get data from GEE and it's working well.

Where I can see Quota usage? There is also no info of using quota in this page. Please see attached file.


texts: Calculate average elevation value in LAS files in a specific area

texts:
I need to calculate the average elevation and its standard deviation of some buildings and parts of roads using LAS files.
I want to be able to consider only the points that are within a polygon in a shapefile. My LAS files have a point space of 0.7 meters and the area covered by each tile is 1500 meters X 1500 meters, but the buildings that I want to consider have small areas of 190 meters X 50 meters.
I only want to have the average elevation of the points that are within the purple box


texts: Calculate a Field in Model Builder in ArcGIS when the field name the calculation is based on is dynamic

texts:
I am essentially trying to do what is listed in this question:

Renaming field via ModelBuilder?

So I am trying to use the Add field > Calculate Field > Delete Field method listed but the problem is my field names keep changing as they are based on an ID (%Value%) in my model so I can't just do

[NewFiled]=[OldField]

Instead I need to do something like

[NewField]= [OldField %Value% +3 ]

However I'm not sure if this is possible to do a calculation within a field name.

Below is the long version:

I have a model (see below) that creates a number of different shapefiles that are named "School_%Value%" where %Value% can be a number from 1 to 797 (I only have 121 shapefiles but their values fall within that range).

Because of the join that occurs half way through the model all the field names get renamed, and because the ID is numeric each field gets iterated so if the Shapefile "School_1" had columns Name, Date, Address, they become School_1, School_2, School_3. If my shapefile "School_100" had the same three columns they will be School_100, School_102, School_103.

I at least know the column I need to change is always the last one, so the name will always be %Value% + 3, but I'm not sure how to do arithmetic in a field value in the calculate field tool.


texts: I have two feature classes. One is a point feature class, the other is a line feature class. The point features all intersect the line feature. I want to route the line feature (have line symbology point to) to point features whose attributes are of lesser value. I havent been able to find an automated option setting up a network, using topology, or using the route tool. Does anyone have any ideas on how to do this?


texts: Scraping Uber service coverage from Uber City website into shapefile?

texts:
I would like to get a full national wide coverage of Uber service area. I find there are specific coverage boundaries by city on its website. For example, https://www.uber.com/cities/charlotte/.

Is there a way extract the polygon and turn it into a shapefile?

If so, is it possible to scale it up for all the cities?


texts: Count Number of Times 100 Feature Classes Overlap

texts:
This must have been asked before, but I can't find the exact answer and my attempts so far have not worked. Most of the overlap questions I see are for calculating the area of overlap for portions of polygons or points in polygons. I want to do something that should be simpler: count overlapping polygons from separate feature classes.

I have 100 separate feature classes, each representing the occurrence of a single species. They all use the same polygon template (based on sub-watershed boundaries), meaning that when more than one species is found in an area the polygon (sub-watershed) for those species' feature classes are perfectly aligned (overlap by 100%, exactly). Note: although not as efficient for disk storage, this is how the data was provided, rather than relating tables of species data.

Example of overlaps: Each color (blue, red dots, orange lines) each represent one species. When they overlap, e.g. red and orange, there is a perfect overlap of polygon boundaries.

What I want to do is count the number of species recorded for any given location, i.e. count the number of times a polygon overlaps and associated this with a specific polygon location. The area of overlap is not relevant in this case, just the count. With a raster, this would be equivalent to setting each layer's cells to have a value of 1, then summing. Is there a simple way to efficiently batch/automate this with 100+ vector features, without converting to raster?

So far I have tried using this count overlaps model.
However, this resulted in an m-aware error, which seems unrelated. I am now running union, which has been progressing for a few hours without completing. Since I may try and run this with 1,000+ feature classes in the future, I was trying to find a more efficient way to complete the analysis. There must be a very simple way to do this, something that I am overlooking.


texts: How to force WMS to return actual image (not blank image) at all zoom levels?

texts:
Some WMS (especially the ones made through ArcGIS) only returns actual images with data at certain zoom levels. At other zoom levels, it will return blank images. I know the server is setup to do that but is there a way around this to force WMS to return images with data at all zoom levels?
The server returns actual image data at this zoom level:

https://hazards.fema.gov/nfhl/services/public/NFHLWMS/MapServer/WmsServer?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=-98.12492909258571672,28.30800648431628375,-98.07080887354277365,28.3716343094073089&CRS=CRS:84&WIDTH=666&HEIGHT=783&LAYERS=4&STYLES=&FORMAT=image/png&DPI=120&MAP_RESOLUTION=120&FORMAT_OPTIONS=dpi:120&TRANSPARENT=TRUE

The same server returns a blank image at higher zoom level:

https://hazards.fema.gov/nfhl/services/public/NFHLWMS/MapServer/WmsServer?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=-98.15361443391026341,28.28411557681083721,-98.04537399582437729,28.41137122699289463&CRS=CRS:84&WIDTH=666&HEIGHT=784&LAYERS=4&STYLES=&FORMAT=image/png&DPI=120&MAP_RESOLUTION=120&FORMAT_OPTIONS=dpi:120&TRANSPARENT=TRUE

(note that it is made using ArcGIS MapServer)
Another server that returns at all zoom level (GeoServer):

https://www.mrlc.gov/geoserver/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=27.76940466477179115,-98.55520921245411614,28.78744986622823632,-97.68928570776702713&CRS=EPSG:4326&WIDTH=666&HEIGHT=784&LAYERS=mrlc_download%3ANLCD_2016_Land_Cover_L48&STYLES=&FORMAT=image/png&DPI=120&MAP_RESOLUTION=120&FORMAT_OPTIONS=dpi:120&TRANSPARENT=TRUE


texts: GRASS: Extend Map Display canvas to full window

texts:
This has been a feature of the 7.8 releases, the map canvas in the Map Display does not extend to the full window.

The size of the map canvas seems to be hard-coded somehow, no matter the size of the window the canvas is always the same size. This is not always critical, but with large maps it becomes difficult to visualise the whole thing.
Is there any way to force the map canvas to extend to the full window size? I work on Ubuntu 20.04.


texts: How to overlay roads layer on top of all layers in OpenLayers?

texts:
I want to put roads layer (from any source, please recommend some..) on top of my KML polygons and base layer:

baseLayer = new OpenLayers.Layer.WMS("OpenLayers WMS", "http://labs.metacarta.com/wms/vmap0", {layers:"basic"});
            map.addLayer(baseLayer);

Is that possible?


texts:
I have the following 2 points on the map: point A lat1 / long1 and point B lat2 / long2 in OL map  I have distance, arrived_time, speed.
With these data how can I get a estimated latitude, longitude  to the next point on 10,20,30 and 40 minutes?
Having start and end point.


texts: WMS on ArcGIS Server 10.6 does not show image tiles but tile footprints

texts:
I have set up a Map Service on ArcGIS Server 10.6 and turned on WMS capability for it.
When I tried to add this WMS to ArcMap 10.6 running on the same machine, it didn't return the actual images but tile footprints.
I also tried add a WMTS service from ArcMap and this worked and showed the image files successfully.
Is it a bug or something I didn't do correctly in setting up the WMS service?


texts: Replace Statistics_analysis with SearchCursor?

texts:
I'm using ArcGIS 10.3, 64bit with an advanced license. I have the following code, that runs just fine. It goes through a table that contains a large number of entries with zip-codes and population information. There are multiple entries for each zip-code. I now want to:

Go through the whole table and select each unique zip code one by one and
calculate the mean value of the population entries

I therefore created the following script:

import arcpy

# set workspace
workspace = r"Path_to_my_WS"
arcpy.env.workspace = workspace
arcpy.env.scratchworkspace = workspace
arcpy.env.overwriteOutput = True

# set inputs
inputTable = r"Path_to_my_table"
field_zip = "Fieldname_with_zipcodes"
field_pop = "Fieldname_with_population"

# create list with unique zip-codes
zipcodes = [row[0] for row in arcpy.da.SearchCursor(inputTable, field_zip)]
uniqueZipcodes = set(zipcodes)

# MAIN TASK: SUM SALES
# empty list to append values to
pop_avg = []
arcpy.MakeFeatureLayer_management(inputTable, "inputTable_Lyr")
# Loop through unique zip codes and perform statistics analysis
for each_zipcode in uniqueZipcodes:
    where_clause = "field_pop =" + each_zipcode
    arcpy.SelectLayerByAttribute_management("inputTable_Lyr","NEW_SELECTION", where_clause)
    statistics = """{0}\statistics""".format(workspace)
    arcpy.Statistics_analysis("inputTable_Lyr",statistics,[[field_pop,"MEAN"]])
    inquiry_pop_avg = arcpy.da.SearchCursor(statistics,"MEAN_field_pop")
    pop_avg.append(inquiry_pop_avg)
    del statistics
print(uniqueZipcodes)
print(pop_avg)

It does work like this, however it takes a rather long time. Would it be possible to replace the arcpy.Statistics_analysis() function with a more performant version of the arcpy.da.SearchCursor() in order to avoid creating and deleting a table for each loop?


texts: Assign lineEdit to Variable for a QGIS 3 plugin

texts:
Good day to all, I am trying to create my first QGIS 3 plugin, I am using Plugin Builder and Qt Designer, I want to add a variable named name with lineEdit to be able to use the information entered later, however, I receive the following error:

name = self.dlg.lineEdit.text()
AttributeError: 'Etiquetas_layout' object has no attribute 'dlg'

here's my code:
def run(self):
    """Run method that performs all the real work"""

    name = self.dlg.lineEdit.text()

    # Create the dialog with elements (after translation) and keep reference
    # Only create GUI ONCE in callback, so that it will only load when the
    plugin is started
    if self.first_start == True:
        self.first_start = False
        self.dlg = Etiquetas_layoutDialog()

    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # Do something useful here - delete the line containing pass and
        # substitute with your code.
        pass


texts: Import OSM with osm2pgsql - missing features

texts:
My problem is missing features in planet_osm_polygon table. In the planet_osm_rel table I can find Sahara desert - it means that relation should exist, but in the planet_osm_polygon table I can't find Sahara at all.  In my planet_osm_polygon table exists hundreds of "relational" features with common osm_is value. For import I use command as below:
osm2pgsql  -H myHostIP  --slim --hstore --create  --number-processes 24 --database osm --username osm  --password --style /openstreetmap-carto/default.style --extra-attributes  --flat-nodes /home/myFolder/flat_nodes/nodes.cache  /home/osm/Pobrane/planet-201109.osm.pbf


texts: bring back column info in sde descendant version

texts: I have a feature class registered as versioned in sde with editor tracking turned on. In my edit version, I was updating one field for all records but forgot to turn off the editor tracking before editing. So the editor tracking fields are updated under my user name and same edit date in my edit version. I have not posted my edit version to the master version so everything in default version is good and there are still new editing going on in default version. Since we really want to preserve the editor tracking info and not to mess up with the edits, what is the best way to fix this problem?


texts: Convert a tiff Float32 image to 8bit tif (or jpeg) image using GDAL and Python

texts:
I need to convert a Float32 tiff image to 8bit format in order to display the image in a web page.
In order to do that I need a jpg image but this conversion is not correctly done with GDAL because the image is Float32 type. So I want to convert to 8bit type with Python.

Can you help me?

I've tried something like this:

gdal.AllRegister()
src = gdal.Open(str(input_file))
src_data = src.ReadAsArray()
cols = src.RasterXSize
rows = scr.RasterYSize
driver = gdal.GetDriverByName("GTiff")
dst = driver.Create(str(output_file), cols, rows, 1, gdal.GDT_Byte)
dst.GetRasterBand(1).WriteArray(src_data)
#dst = drv.CreateCopy(output_file, src, 0, saveOptions)

I try with Create() and CreateCopy() but I always get a black image.


texts: Store Encrypted DB credentials in QGIS Server

texts:
I have deployed a QGIS server following steps in the user manual. https://docs.qgis.org/3.16/en/docs/server_manual/containerized_deployment.html#docker-stacks
I am serving some vector layers as WFS using the QGIS project on the server. The vector layers are stored in a PostGIS DB. While working on a desktop one can save the credentials encrypted in 'qgis-auth.db'. How can one do the same if the project is on a server.
The server is running on a Linux EC-2 instance.


texts: Filtering image collection based on list of dates using Google Earth Engine Python API

texts:
I would like to extract an image collection of some days.
For instance, if I want to extract an image collection between two dates, I can simply do:
START = '2010-08-01'
END = '2010-09-15'
cover = ee.ImageCollection('MODIS/006/MYD11A1').filter(ee.Filter.date(START, END))

What I would like to do is to filter based on a list of dates and do something like the following:
DATES = ['2010-08-01', '2010-08-12', '2010-09-01']
cover = ee.ImageCollection('MODIS/006/MYD11A1').filter(ee.Filter.date(DATES))

Is it possible?


texts: How to change OpenStreetMap Layer Properties in QGIS 2.0

texts: I have loaded the OpenStreetMap layer in QGIS 2.0, but would like to edit the properties in order to better complement the map I'm producing. I can't seem to find the OSM toolbar or OSM Feature widget that allows me to edit. Is this  because OSM has been integrated into the core functionality of QGIS 2.0?


texts: How do I get the latest version of QGIS to find my plugins

texts: I just downloaded QGIS 3.12.2-București for the Mac (MacOS Mojave 10.14.6, 3.2 GHz with 32 GB ram).  I've been running an older version (QGIS 2.18.2-Las Palmas) for a couple of years and have quite a few plugins for it.  How do I get the new version to locate and recognize the plugins I have already downloaded?


texts: Setting transparency for no data in raster in geoserver

texts:
I have a raster published in my geoserver with  default raster style as below

<NamedLayer>
<Name>default_raster</Name>
<UserStyle>
<!-- Styles can have names, titles and abstracts -->
  <Title>Default Raster</Title>
  <Abstract>A sample style that draws a raster, good for displaying imagery</Abstract>
  <!-- FeatureTypeStyles describe how to render different features -->
  <!-- A FeatureTypeStyle for rendering rasters -->
  <FeatureTypeStyle>
    <Rule>
      <Name>rule1</Name>
      <Title>Opaque Raster</Title>
      <Abstract>A raster with 100% opacity</Abstract>
      <RasterSymbolizer>
        <Opacity>1.0</Opacity>
      </RasterSymbolizer>
    </Rule>
  </FeatureTypeStyle>
</UserStyle>

And my imagery looks like

Which has black color set for no data. I want to make the no data to transparent. How that can be done?


texts:
I'm trying to automate a map, and one of those tasks is to set a specific symbology for a feature class. I create the feature class of Subject and Comps in the code, then add it to the map, then when I try to add the symbology created previously, it gives me an error

ERROR 000732: Input Layer: Dataset C:\ArcMap\Res_Parcels.gdb\Subject does not exist or is not supported
Failed to execute (ApplySymbologyFromLayer)

arcpy.env.overwriteOutput = True
# Set layer to apply symbology to
inputLayer1 = "C:\ArcMap\Res_Boulder.gdb\Subject"
inputLayer2 = "C:\ArcMap\Res_Boulder.gdb\Comps"
# Set layer that output symbology will be based on
symbologyLayer1 = "C:\ArcMap\Comps.lyr"
symbologyLayer2 = "C:\ArcMap\Subject.lyr"
# Apply the symbology from the symbology layer to the input layer
arcpy.ApplySymbologyFromLayer_management (inputLayer1, symbologyLayer1)
arcpy.ApplySymbologyFromLayer_management (inputLayer2, symbologyLayer2)


texts: pgRouting: Find source/target ID of edge closest to points in a separate layer for many-to-one routing

texts:
I have a small road network which I have extracted from OSM and prepared for routing using OSM2PO. I also have a separate layer of points (postcode centroids). I wish to compute many-to-one routing between each point to a single destination.

In my research I have found many useful answers regarding how to undertake a many-to-one routing, but all these solutions require manually entering a source/target id for each origin edge. (The single destination is, of course, simple and I am happy to manually input that source/target id.)

Is there a way in which I could programatically find the nearest edge to each point and then pass the relevant source/target ID into the pgRouting algorithm, looping through for every point?

Here is a picture of my data in QGIS which may help clarify the situation:


texts: In the UK, how do I convert from an address to a lat/long?

texts: Please include where I can buy the dataset and how much it costs?


texts: gp.CreateFeatureClass(..) and gp.CopyFeatures(..) isn't copying all fields of feature class

texts:
How do I get the ArcGIS 9.3 geoprocessor to copy all of the fields of a feature class?

Here is a description of what I'm trying to do:

I am modifying an existing script created by a developer not on the program anymore. The script's purpose is to extract a given feature set into a zipped .gdb.  Therefore, the script takes a feature set as input. In the FeatureSet that I am attempting to export (which is a feature set representing major lakes), the fields are DEPTH, NAME, OBJECTID, Shape_Area, Shape_Length, SQKM, SQMI, and SURF_ELEV. The .gdb resulting from this script, however, only contains OBJECTID, Shape (where'd that come from??), Shape_Length, and Shape_Area.

I am exporting the feature set via

gp.CopyFeatures(featureSet, scratchGdb + os.path.sep + lyrName)

The .gdb is successfully created, but only some fields are there!  Why am I losing some of the fields in this feature class when exporting it to a .gdb?

Through using Komodo IDE, I was able to remote debug the script, and I saw that the input feature set only had the four fields that ended up in the exported .gdb.

Does this have to do with setting a schema for the feature set?


texts: Display full date-time, even when the time is 00:00:00

texts:
I have date fields in Oracle tables. When I view the data in the attribute window, the dates are displayed like this:

1/1/2017
1/2/2017 3:50:28 PM

As demonstrated in #1, ArcMap doesn't show the time portion of the date if it is 00:00:00 (even though it technically exists in the underlying table).

Is there a way to force ArcMap to show the full date-time of date fields in the attribute table, even if the time component is 00:00:00?

Example:

1/1/2017 00:00:00  <<<---
1/2/2017 3:50:28 PM

I've poked around the documentation, but haven't found anything specific to this exact issue.

Fundamentals of date fields
Changing the way ArcMap displays short format dates in Windows Vista

The reason I ask is:

I'm doing some intensive testing on views based on temporal tables. I'm doing a lot of querying and editing of the time portion of date fields. It would be helpful to see the time portion of the dates, even when the time is 00:00:00.


texts:
I'm running ArcMap 10.3 and am wondering the quickest way to export a map from the layout view without navigating up to the File menu & clicking the "Export Map" feature. Which is to say, I'm looking for a keyboard shortcut for this feature.

I can't find a built-in keyboard shortcut for exporting but I've noticed that if the File drop down menu is clicked, I can simply hit the E key and that will bring up the Export Map dialog box. Does anyone know of a shortcut to at least open/expose the File drop down menu in the top left? If so, then it's a quick 2 strokes and I'm there.


texts: self.iface.mapCanvas().destinationCrs() Replacement in QGIS 3

texts:
I need to store the map canvas' destinationCrs in a variable, but when using the following line:

self.iface.mapCanvas().destinationCrs()

I get the error:

  'QgsMapCanvas' object has no attribute 'destinationCrs'

What is the appropriate replacement for this function in QGIS 3.x?

SOLUTION:  The REAL issue was omitting the MapRenderer() instead of replacing it with mapSettings(), which appears to have resolved the dilemma.

Original (QGIS 2) code:

    map_crs = self.iface.mapCanvas().mapRenderer().destinationCrs()

New (QGIS 3) code:

    map_crs = self.iface.mapCanvas().mapSettings().destinationCrs()


texts: How to nest mapped functions with the Earth Engine Python API

texts:
I'm trying to extract data from multiple Images in an ImageCollection using multiple Features in a FeatureCollection. Because the Features are spread over a large geographical area, I can't supply the whole FeatureCollection to the reduceRegions() function without running out of memory. Instead it is easier to map over the FeatureCollection, and map the reduceRegions() function over the ImageCollection. Using flatten() at each map then results in a nice table for export to CSV. In JavaScript we can use anonymous functions, so it is relatively straightforward  to map a function within a mapped function as such:
var img1 = ee.Image(1);
var img2 = ee.Image(2);
var imgCol = ee.ImageCollection.fromImages([img1, img2]);

function mapReducerOverImgCol(feat) {
  function imgReducer(img) {
    return img.reduceRegions({
      reducer: ee.Reducer.first(),
      collection: ee.FeatureCollection(feat),
      scale: 500,
      tileScale: 10
    });
  }

  return imgCol.map(imgReducer).flatten();
}

var table = featCol.map(mapReducerOverImgCol).flatten();

print(table);

https://code.earthengine.google.com/b031ed7d8037181a4398c38b00dd16de
However, I'm trying to figure out how to do this in Python. Here is my failed attempt:
def nestedMappedReducer(featCol, imgCol):

    def imgReducer(img):
        return img.reduceRegions(
        collection = ee.FeatureCollection(feat),
        reducer = ee.Reducer.first(),
        scale = 500,
        tileScale = 10
        )

    def mapReducerOverImgCol(imgCol):
        return imgCol.map(imgReducer).flatten()

    return featCol.map(mapReducerOverImgCol).flatten()

However this fails and complains that 'feat' is not defined.
Any ideas on how to get this to work?


texts: Adding file geodatabase attachments to features using QGIS

texts:
Is it possible to add attachments (PDFs, Word documents, etc) to features in QGIS (and not via a hyperlink)?
An attachment, as in a binary file that is stored in a file geodatabase, rather than a hyperlink.


texts: Getting access to field using its name in PyQGIS

texts:
I have two layers: Data and Thresholds. Data is a point vector layer (GeoPackage) that has hundreds of fields. Thresholds is a csv table with three threshold values for each field. My goal is to create a new field in Data that would be a classification value. I need to check if any value of a feature in Data exceeds the threshold value of the same field in Thresholds. For example, if any Data value exceeds the first threshold value in the corresponding field in Thresholds, the classification value would be 1. Or if it exceeded the third threshold value, the classification would be 3. So I will have to loop through the values in Thresholds for each field in each object in Data. A couple nested loops, no biggie here.

The way I want to do this is to take the name of the field in Data and check if there is such a field in Thresholds. If there is, then I will loop through the field and change the value in Data. It is possible that in the future, Data will contain fields that are not in Thresholds.

How do I get access to a fields name, and how do I get access to a field using just the name?


texts: Getting layer id of processing step in QGIS graphical modeler for usage within aggregate() function

texts:
How can I get the layer name or layer id of the result of an processing algorithm within QGIS graphical modeler, so I can use aggregate() function in a later processing step? Is this somehow possible at all?

aggregate( @layer_name,'collect',$geometry) cannot work here because @layer_name is not related to a specific processing algorithm output. If I try aggregate( @name_of_processing_step_OUTPUT,'collect',$geometry) or aggregate(layer_property( @name_of_processing_step_OUTPUT,'id'),'collect',$geometry) I get Cannot find layer with name or ID '' in preview and just no output when running the model.


texts:
I want to use the following Leaflet function:

# Generated using folium

<!DOCTYPE html>

    <meta http-equiv="content-type" content="text/html; charset=UTF-8"/>

    <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/leaflet@1.5.1/dist/leaflet.css"/>
    <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.2.0/css/bootstrap.min.css"/>
    <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.2.0/css/bootstrap-theme.min.css"/>
    <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/font-awesome/4.6.3/css/font-awesome.min.css"/>
    <link rel="stylesheet"
          href="https://cdnjs.cloudflare.com/ajax/libs/Leaflet.awesome-markers/2.0.2/leaflet.awesome-markers.css"/>
    <link rel="stylesheet"
          href="https://rawcdn.githack.com/python-visualization/folium/master/folium/templates/leaflet.awesome.rotate.css"/>

    <meta name="viewport" content="width=device-width,
                initial-scale=1.0, maximum-scale=1.0, user-scalable=no"/>

    <link rel="stylesheet"
          href="https://cdnjs.cloudflare.com/ajax/libs/leaflet-locatecontrol/0.66.2/L.Control.Locate.min.css"/>

I want that the location found when this function is activated will be saved, maybe with AJAX or to JSON.

JSON file would be ideal since I want to work with the data in python as well.

I'm kinda new to JavaScript so my knowledge on this click on works is limited (I mostly work with Python).


texts: Create random points within polygons of a spatial data frame in R

texts:
I want to create random points within polygons that are stored in a spatial data frame in R.

I tried to use the csr function of the splancs package, and the spsample of the sp package without success. I typed the following code:

Csr(Greece.map[1, ], 12)

And received:

  Error in storage.mode(poly) <- "double" no method for coercing this s4
  class to a vector

And sample(Greece.map[1, ], n=12, random) receiving the following error:

  Error in sample.Spatial(as(x, "Spatial"), n_tot * (1 + its * 0.1),
  type = type, : object random not found.


texts: Calculating azimuth of direction in OpenLayers?

texts:
I am developing an application with OpenLayers and I need to calculate the angle between two segments or directions.

I thought about getting the azimuth of those segments and then calculate the difference but I cannot find any method that can perform this directly, something like ol/geom/LineString getInitialAzimuth() or similar.

I found this http://www.movable-type.co.uk/scripts/latlong.html that could be useful but I was expecting something simpler.

I am missing anything?


texts: QGIS 2.6 how to prevent user to drop a postgis table

texts:
Does anyone know how to prevent user to drop a postgis table?
I have a vector layer imported from a postgis2.0 database, and one of my users drop it accidently. Fortunate I had a backup, but i don't expect that to happen again.
we're working on Win7 64bits, Postgresql 9.4 - Postgis2.0 and Qgis 2.6.
best regards.


texts: How to change the resolution of maps exported from QGIS Time Manager

texts:
Is it possible to change the default resolution of the 'Save as image' export option - which is used by the Time Manager plugin - to export higher quality images that can be used in animations/movies?

I've dug deep into the settings and even the ones I'm warned not to go into, but can't find anything.


texts: Correcting commercial grade GPS data to 2 meter accuracy in RT?

texts: I have a set of rovers with commercial grade GPS boards (Adafruit Ultimate GPS Breakout). They send their uncorrected data to the central server. There, before recording each position in the database, I need to correct it the best I can, hopefully to a 1-2 meter accuracy. I'm new to this topic and have read a lot of info on sending corrections to the rover's GPS units and have them correct itself. I don't need this, as I'm correcting everyone in a central server location. I imagine, this should be just a matter of adding a small number to every long, lat, and elev reading I've got from the field. But from what I've read so far, no one is doing just this. Is there a simple way to get these corrections and apply them in this simple way? I can set up a base station (another Adafruit Breakout?) that would know it's location and would send it's reading to the server as well. Would it be enough to just get its reading, compare to it's known location, and apply the 3 deltas to corresponding readings from the rovers in the field?


texts: How to use 'remap' in reclassify tool code in ArcGIS Python Toolbox

texts:
I am making a reclassify tool in ArcGIS python toolbox (.pyt) ...
I am able to use the remap datatype but am not able to parse the information present in Remap in the Reclassify function's third argument ...

Code is below:

import arcpy, string
import numpy
from arcpy import env
from arcpy.sa import *

class Toolbox(object):
def __init__(self):
    """Define the toolbox (the name of the toolbox is the name of the
    .pyt file)."""
    self.label = "Toolbox"
    self.alias = ""

    # List of tool classes associated with this toolbox
    self.tools = [Remap]

class Remap(object):
def __init__(self):
    self.label = "Remap"
    self.description = ""
    self.canRunInBackground = False

def getParameterInfo(self):
    feature_class = arcpy.Parameter( displayName="Input Raster ", name="input1", datatype="Raster Layer", parameterType="Required", direction="Input")
    field = arcpy.Parameter( displayName="Input Field ", name="input2", datatype="Field", parameterType="Required", direction="Input")
    remap = arcpy.Parameter( displayName="ReMap ", name="remap", datatype="remap", parameterType="Required", direction="Input")
    out_class = arcpy.Parameter( displayName="Output Raster ", name="ouput1", datatype="Raster Layer", parameterType="Required", direction="Output")

    field.value = "Value"
    field.enabled = True

    remap.parameterDependencies = [feature_class.name,field.name]

    params = [feature_class,field,remap,out_class]
    return params

def isLicensed(self):
    return True
def updateParameters(self, parameters):

    #if parameters[0].value:
        #parameters[1].parameterDependencies = [parameters[0].name]
        #parameters[2].parameterDependencies = [parameters[0].name,parameters[1].name]
    return

def updateMessages(self, parameters):
    return

def execute(self, parameters, messages):
    inRaster = parameters[0]
    outRaster = parameters[3]
    remap2 = parameters[2]
    arcpy.gp.Reclassify_sa(inRaster,"Value",remap2,outRaster,"DATA")
    return


texts: How to use coding for auto-comparing two different fields in ArcGIS for Desktop?

texts:
I want to compare two group fields in attribute table, which are related to the name of villages.

The name of these villages (in both field) are the same but they have in some cases differentiation (because of misspelling).

In order to update my databases I need to find different names and edit them.

They are many names and it doesn't make sense to check them manually.

I think Python or VBA code in ArcGIS is practical.

How can this be done?


texts: ArcGIS Calculate Value tool - Python syntax error

texts:
I have a feature class within a file geodatabase.

I am using the "Iterate row selection" iterator in Model Builder to read each row in the FC.
I am then using the "Get Field Value" tool to read the following 3 fields from each row:

"Status"

"User"

"Tile"

I then want to use the "Calculate Value" tool to return a string, which I will later use as a file path to copy a specific file.
The logic is as follows:
If "Status" is "Available", return C:\Data\ + "Tile"
otherwise return C:\Data\ + "Tile" + _ + "User"

In the Calculate Value tool, I've put the following into the expression field:

a("%Status%", "%User%", "%Tile%")

I've put the following into the Code Block:

def a(Status, User, Tile):
    if Status == "Available":
        return "C:\Data\" + Tile
    else:
        return "C:\Data\" + "Tile" + "_" + User

I have the "Data type" set as string.

My model is as follows:

I get the following error for my code block:

    ERROR 000989
Python syntax error: Parsing error SyntaxError: EOL while scanning string literal (line 3)

I'm happy to have a solution using model builder or using a stand along python script.


texts: Converting Lat Lon in WGS84 format to UK NGR in QGIS

texts:
I have csv file that has WGS Lat/Lon, what I need to do is convert these to UK NGR i.e Easting/Northing.
In QGIS there under Field Calculator I think there are some defined operator like $X and $Y that does this but I can't make it work.


texts: How to use ArcGIS 10.2 Features to JSON output in web maps?

texts: Using the Features to JSON tool that is now available with 10.2, is there additional formatting that needs to be done to the files before they can be displayed as web maps?


texts: How to create an if statement for moving data if a cell is blank from another column

texts: I need help creating an if statement for moving data from one column to another if it is blank. Any suggestions? See screen shot for highlights noting the data that needs to be moved (po boxes).


texts: Can anybody tell me why there is still a python27 folder along with the python37 folder in the latest (3.x) QGIS releases?


texts: Raster Calculator not excluding 0 values (QGIS)

texts:
I am using QGIS 3.10.2 and attempting to mask out 0 values with the Raster Calculator. I have a series of values from 0 to 255, with 0 being empty space on the boundaries that I do not care about. My equation in the Raster Calculator is as follows:

("a">0)*"a"

Where "a" is the raster layer I am trying to mask. However, the output gives me the same raster, with 0 values still included. Can anyone shed some insight on this problem?

I am new to QGIS.


texts: How to fix positive North-South Pixel size in GeoTIFF with GDAL

texts:
I have issues with some GeoTIFF raster files with positive NS resolution (North-South)

gdalinfo gave me this:

Pixel Size = (1.000000000000000,1.000000000000000)

I run into all kind of problem with gdal_merge or gdalbuildvrt when trying to mosaic a set of files.

I get error messages like:

gdalbuildvrt does not support positive NS resolution

The solution should be to flip the whole thing by moving the vertical origin to other side of the N-S extent and invert Pixel Size/resolution to a negative value, as almost all GDAL tools expect:

Pixel Size = (1.000000000000000,-1.000000000000000)

How is this best achieved? (Preferably with exiting GDAL command line tools)


texts: What's the correct way to reproject a population raster?

texts:
I'm working with a population raster that is in the WGS84 coordinate system, and I want to get it into a projected coordinate system (Albers Equal Area), using ArcGIS 10.1.

According to the metadata, the grid cells are "30 arcseconds (1 km)" in size (on a side presumably) and their values represent the population count in that 1 km^2 area.

When I use Project Raster, the tool automatically fills the "Cell Size" parameter with a value around 830 meters. When I compare equivalent 3-by-3 windows of data cells between the WGS84 and the Albers versions, I'm finding that values have been changed.

Here's WGS84 (ignore the 4th column of cells for now):

Here's the Albers equivalent of the 3-by-3 window above:

Note that, when the WGS84 one is projected to Albers, the values in WGS84's 3rd column are "displaced" by the values in the 4th column. I'm assuming it's doing some kind of nearest-neighbor resampling because of the altered cell size? The ranges on the 2 rasters also differ by ~2,000 on the top end.

Further points of confusion:

How can the metadata claim that 30 arcseconds = 1 km when the distance of an arcsecond varies depending on longitude? Is it referring to the conversion at the equator? I was under the impression that at the equator, 1 arcsecond = 30 m, so 30 arcseconds = 900 m. Is the metadata simultaneously implicitly stating "at the equator" AND rounding up?
In ArcGIS's Project Raster tool, where is the default cell size of 830.4934540 coming from? Does this mean each cell is to be 830 m on one side, which would be about 0.69 km^2? Should I alter this to be more in keeping with the metadata's value of 1 km^2 per cell?

What is the correct way to do this reprojection while remaining loyal to the original population counts? Why are my values clearly being altered so there is no one-to-one correspondence and what's the best practice for accounting for this?


texts: Seeking quarternary watershed map of Eastern Ontario?

texts:
I am looking for a quarternary watershed map of Eastern Ontario.

Where is the best place to look for this?

I am currently using the ArcGIS Online basemap software to create my maps. I have a CSV file with the watersheds of Eastern Ontario and their geographical polygon boundaries but arcGIS online cannot read it.


texts: Problem reprojecting LIDAR data with Liblas

texts:
I am currently trying to reproject a LIDAR file (.las) using the Open-Source LibLAS, compiled on a Scientific Linux server.
The Original data is in Oregon Lambert, NAD83 (EPSG:2992), and I am trying to reproject it to UTM WGS 84, Zone 10 N (EPSG:32610).
However, when I use the --t_srs option, I get the following error:
"error: X scale and offset combination is insufficient to represent the data"
I am suspecting it has to do with the integer values stored and the distance of translation.
How do I get it to work? Is there anything I have to change in the header information?
Regards,
Edit: Here is the command line I am using (lasinfo detects the original EPSG, so no need to specify it):

./las2las -i path/to/input/folder/*.las -olas -odir path/to/output -v
--t_srs EPSG:32610

Result:

$ Setting output SRS to EPSG:32610
$ error: X scale and offset
combination is insufficient to represent the data


texts: Labeling from a defined string in MapBasic

texts:
I am currently working on a code to apply a variety of different changes to layers with certain attributes, I have been successful with this so far (thanks to Peter). However I now want to add a text box to this so the user can enter the name of the field they want labeling (as this is the only thing differs in the .TABs I intend to use this tool for).

I have currently created a dialog box for the user to input this however when it comes to labeling (I assume), MapInfo is looking for a field with the same name as my string to label with, instead of looking at what is stored there.

For example

Dim sLabinput As String

Set Map Window nMID
            Layer nLayer
                Display Global Global Line (2,2,16711680)
                Arrows On
                    Label With sLabinput Auto On Visibility On
                    Font ("Arial",257,9,0,16777215)

I haven't put the rest of the code as I want to try and learn this myself however if that is the only way to solve it I can post it.

Am I right in assuming in this instance MapInfo is looking for a field called sLabinput, rather than the value stored in there?

Thanks,
Alan


texts: NameError: name data_all1 is not defined in Python 3.x

texts:
I have a standalone script using data opened from a Geopackage with Fiona. When I run this portion of my script:

import os
import fiona
import numpy as np

datapath= '/media/sf_SharesShapefiles/'

layer_keys =fiona.listlayers(os.path.join(datapath, "test.gpkg"))

for l in layer_keys:

    layer1 = fiona.open(os.path.join(datapath, "test.gpkg"),layer=l)

seasons = ['ann','djf','mam','jja','son'] RCPS = ['rcp45', 'rcp85'] vari =  ['Precipitations_totales_CRU_3_21TS_CMIP5_']

nColors = 10 # number of categories Eq_break_flag = False # equal breaks vs quantiles NoData =-9999 # nodata value if required

break_vals_Seas ={} for v in vari:
    break_vals_Seas[v] = {}
    for s in seasons: # loop season\
        for r in RCPS:  # loop rcps

            layer1 = fiona.open('/media/sf_SharesShapefiles/test.gpkg', 'r', layer = 'Precipitations_totales_CRU_3_21TS_CMIP5_' + r)
            #column_names = [row[:] for row in cols]
            dataorig = {}

        for i in range(0, len(layer1)):

            rec = layer1.next()
            test = rec['properties']

            if i == 0:
                for k in test.keys():
                    dataorig[k] = np.asarray(test[k])
            else:
                for k in test.keys():
                    dataorig[k] = np.append(dataorig[k], test[k])

        # print(data1.keys())
        column_names = list(dataorig.keys())

        indices = [i for i, s in enumerate(column_names) if 'CRU' in s]
        #df = (dataframe[:, indices[0]:3])
        column_names = column_names[indices[0]:]

        for c in column_names:
            if 'Shape' and 'Length' in c:
                column_names.remove(c)

        fieldNamesSubset = [ss for ss in column_names if s in ss]
        for ff in range(0,len(fieldNamesSubset)): # Concatenate data from field subset (each season) into a single numpy array

            data1 = dataorig[fieldNamesSubset[ff]]  #arcpy.da.FeatureClassToNumPyArray(l,fieldNamesSubset[ff]) # read field into numpy array

            data1[data1==NoData]=np.nan
            if ff == 0 and r in RCPS[0]:
                data_all1 = data1
            else:
                data_all1 = np.vstack((data_all1,data1))

    if Eq_break_flag:
        interv1 = (np.max(data_all1[~np.isnan(data_all1)]) - np.min(data_all1[~np.isnan(data_all1)])) / nColors
        break_vals = []
        break_vals.append(np.min(data_all1[~np.isnan(data_all1)]))
        for i in range(1, nColors):
            break_vals.append(break_vals[i - 1] + interv1)

        break_vals.append(np.max(data_all1[~np.isnan(data_all1)]) + .001)

    else:
        interv = 100 / nColors
        np.arange(0.0, 1.0, 0.1)
        break_vals = np.percentile(data_all1[~np.isnan(data_all1)], np.arange(0, 100 + 1, interv))

        break_vals_Seas[v][s] = break_vals
        print(break_vals)

I get this error message:

break_vals = np.percentile(data_all1[~np.isnan(data_all1)], np.arange(0, 100 + 1, interv))

NameError: name 'data_all1' is not defined

I think this might be a scope error, however I am not sure how and where to define data_all1


texts: QGIS 2.0 portable or silent installation

texts:
Is it possible to install QGIS 2.0 unattended (silent) on Windows-Client (XP,7) via commandline ?

Or does there exists a portable version for 2.0 ?


texts: Free GIS software for mapping drain lines on a golf course

texts:
I work on a golf course that is built on old farm land. We have old clay drains from the turn of the 20th century and lots of additional drains that have been added in the meantime years.

No maps have been produced and we are constantly discovering forgotten work. Our problem is that we cannot remember all of the lines on which they run so would like to start mapping as we are probing, rodding and digging.

Please can someone suggest a free software that we can use to produce maps?

We just need to plot lines on either maps (or preferable aerial images) of the land with some sort of way of locating them in the future. Either using datum points (trees etc) or GPS (though we will only have our mobile phones to get our co-ordinates.

The golf course in in Wales, UK.


texts: Setting LAS Class Codes in ArcGIS Pro too slow

texts:
I have seen similar questions, but they don't quite relate to my particular issue.
I classified habitat in a high-resolution drone orthomosaic image of a 300-m long stretch of beach, then used Raster to Polygon to convert the classified raster to vector format.  I only have a few different classes: sand and several classes of seaweed.  However, because the seaweed is scattered across the beach in many different shapes and sizes, many many polygons were generated for each class from Raster to Polygon (I used this tool with and without checking the Multipart Feature option, both have the same result).  This has proven to be problematic for my next step: classifying the LAS point cloud.
I have an unclassified LAS point cloud generated from Drone Deploy that I now need to classify according to those habitat polygons.  The Set LAS Class Codes Using Features tool is just the tool I need and I know it works - I just classified my LAS point cloud using one large survey area polygon as the input 'Feature' to classify to and it processed in seconds.
However, this tool is not yet able to finish processing the LAS points using the habitat class polygons, as each class is comprised of many many multipart polygons. By struggling, I mean it has yet to process anything past 13%, even after an overnight run for over 9 hrs. I have tried this tool on ArcGIS Pro 2.5, 2.6, and Desktop, as well as the Extract LAS tool, all with the same result.  The original LAS point cloud was 15 million points, so I also just used the Thin LAS tool to limit the points to ~5 million.  This is currently running overnight on the habitat class with the smallest coverage, so stay tuned.
Is there any way around this problem in ArcGIS Pro?
This is time-sensitive and I have multiple habitat classes to reclassify. Again, this works great one one simple polygon, but not in my case where each class is comprised of thousands of different shaped polygons.  At this rate it will take days to classify the LAS point cloud for one image, time I don't really have.


texts: Why are some GRASS raster tools in QGIS changing the CRS to from EPSG:31287 to USER:100029 for the processed result files

texts:
I am processing raster data for Austria in EPSG:31287 (MGI /Austria Lambert) and all input data as well as the project itself have this CRS properly defined.
When I run different raster tools from the GRASS toolbox in QGIS (I used Version 3.12 and 3.14) (in my case r.reclass and r.neighbors) the CRS of the output file gets a user defined CRS named USER:100029 plus a lot of projection parameters (those parameters are ident to MGI /Austria Lambert specifications plus WGS84 to MGI transformation parameters). The raster cells of this layers are visually slightly shifted in the map view (see images below, the raster cells are 10x10m).
When I change the CRS in the layer properties > source to EPSG:31287 or assign EPSG:31287 as projection (with GDAL > Raster projections > Assign projection) everything is fine. The thing I want to understand this issue and search a for a solution to avoid the additional this projection assignment.
Furthermore, I am teaching GIS, and the rookie students a really confused with problems like that...
Is this a bug, or can this be fixed somehow?


texts:
Is there a way to zoom in a map on an R spatial polygons object.

i.e. if I have polygons for a whole country but I want to zoom in on a specific city and just plot that area how do I do that ?

At the moment I am just plotting a chloropleth of the whole country like this:

    colours<-(brewer.pal(9, "YlGn"))
    brks<-classIntervals(poly$var1, n=9, style="quantile")
    brks<-brks$brks
    plot(ED, col=colours[findInterval(poly$var1, brks,all.inside=TRUE)], axes=F, border=FALSE)

I wanna zoom in and plot a rectangular area around a particular city but can't figure out how.


texts: Updating a feature class based on evaluating different fields with arcpy.da.UpdateCursor

texts:
I have a feature class with a number of fields, two of wich are key fields, named key and front, and secondary fields named key_n and front_n.

The thing is I need to search in the whole table for any instances of key values where the front equals the front_n value (key+front = key+front_n), as the basis for calculating the key_n value.

I have tried with arcpy.da.UpdateCursor, but haven't got it, since, to my view, the arcpy.da.UpdateCursor compares only values within the current row, and the provided examples with ArcGis don't explore this possibility.

This is the code I've gotten lately, to no success. In this case, the fields cvegeo and cveft are to be searched as a match to cvegeo + cveref1, in which case, the cvevial1 and cveseg1 fields are to be carried from the record where cvegeo+cveft are taken originaly.

import arcpy
import sys

frentes = arcpy.GetParameterAsText(0)
#frentes = arcpy.GetParameterAsText(1)
campos = ('CVEGEO', 'CVEFT', 'CVEVIAL', 'CVESEG', 'CVEGEO_N', 'CVEFT_N', 'CVEVIAL_N', 'CVESEG_N', 'CVEREF1', 'CVEVIAL1', 'CVESEG1', 'CVEREF2', 'CVEVIAL2', 'CVESEG2', 'CVEREF1_N', 'CVEVIAL1_N', 'CVESEG1_N', 'CVEREF2_N', 'CVEVIAL2_N', 'CVESEG2_N', 'PRUEBA')

with arcpy.da.SearchCursor(frentes, campos) as scursor:
    for row in scursor:
        abuscar = row[0]
        frente = row[1]
        clave_vial = row[2]
        clave_seg = [3]
        with arcpy.da.UpdateCursor(frentes, campos) as ucursor:
           for urow in ucursor:
               if (urow[0] == abuscar and urow[8] == frente):
                   urow[9] = clave_vial
                   urow[10] = clave_seg
                   ucursor.updateRow(urow)


texts: Spatial interpolation from categorical data in R

texts:
I have a spatial point data set, where each point represents a person pronouncing a word in a specific fashion. Let's say I have six different pronunciations.

In an earlier iteration, I aggregated the point data to a hexagonal grid, the color represents the dominant pronunciation in the hexagon while the opacity represents its dominance, see: https://www.srf.ch/static/srf-data/data/2016/dialektkarten-karte/

Now I'd like to spatially interpolate the data to get a continuous surface, similar as Josh Katz did: http://www4.ncsu.edu/~jakatz2/files/dialectposter.png

He uses something he calls "k-nearest neighbor kernel smoothing" but he's not publishing his R code so I don't fully understand how that is supposed to work with categorical data. I perfectly understand on how to do it with continuous variables, but not with the latter.

The pseudo-algorithm I came up with so far:

Repeat the "final estimate" formula for each pronunciation, using yi = 1 for points where that pronunciation is used, yi = 0 if not. By doing so, you get a grid for each pronunciation. The cells represent the probabilities of the respective pronunciation appearing at that grid cell.
Now, find the most probable pronunciation in each grid cell. This determines the cell color (hue). Its probability determines the color value (i.e.: the higher, the darker, and vice versa).
Find a way to map that resulting grid / raster - might be tricky because it contains two variables that need to be mapped to visual variables (dominant pronunciation and its probability).

But maybe this is completely off - and also, I wonder whether there are R packages for spatial interpolation of categorical data. So far I have only found knn methods and the like and I suppose I have to build something myself.


texts: Determining search variable when using Inverse Distance Weighted Interpolation?

texts:
When using IDW interpolation in ArcGIS, how should I determine the number of nearest input sample points to be used to perform interpolation?

The default number is 12, but it may be incorrect.


texts: Find Point Along Line at Certain Distance

texts:
We have software that returns '5143m from PointA on PathB' and I need to know where that would be without using the measuring tool and snapping to manually trace the line. I am not looking to measure distances between points on a line nor find the distance from points to a line. Is there a plugin I can use or a native feature?

I need to figure the physical location of a point on a line (Point X, Line A) given the point at the start of the line (Point A) and a distance along the line. I.e., if I am given '4000m from Point A on Path B', then manually tracing Path B for 4000 meters from Point A should give me the location of Point X; whatever solution we find here should return the same result, saving me the tedium and time (I have quite a few of these to do).

I've done some more searching via linked questions, and solving this problem in ArcGIS would involve Linear Referencing and Interpolation. Is there something similar for QGIS?

The line I need to be traced is not straight; it is a polyline with vertices, turns, etc.


texts: How to blend Landsat scenes in image mosaic?

texts:
I have Landsat scenes in the same month but different dates.  All are atmospherically corrected. I want to mosaic them but the colors/tones of the final result did not blend. I do not think adjusting the brightness, saturation and contrast would work, as suggested in this post.
I tried QGIS by following GDAL->Raster miscellaneous->Merge but it didn't work. The same result with ArcGIS, the colors didn't blend. Is there any other tool in QGis or ArcGIS to properly do this?
Here's a sample of my result.


texts: Convert CSV file in shapefile with pyshp/shapefile

texts:
I'm trying to convert a file.csv to a file.shp through the pyshp / shapefile module.
The content of the CSV is:

Lon,Lat,Vel,Dir
15.9576198620869825,41.1791119428255143,9.3949995040893555,92.0000000000000000
15.9598119351079415,41.1791119428255143,9.3949995040893555,92.0000000000000000
15.9620040081288987,41.1791119428255143,..., ...

the code I wrote is as follows:

import csv
import shapefile as shp

points = shp.Writer(shp.POINT)
with open(csv_path + '.csv', 'rb') as csvfile:
     csvreader = csv.DictReader(csvfile)
     header = csvreader.fieldnames
     [points.field(field) for field in header]
     for row in csvreader:
        points.point((float(row['Lon'])),(float(row['Lat'])))
        points.record(*tuple([row[f] for f in header]))
points.save(r'path/to/output/file.shp')

I always used this code without problems, but I don't understand why in this case it return:

Traceback (most recent call last):
  File "my_script.py", line 110, in <module>
    points = shp.Writer(shp.POINT)
  File "/home/username/.local/lib/python2.7/site-packages/shapefile.py", line 1018, in __init__
    self.shp = self.__getFileObj(os.path.splitext(target)[0] + '.shp')
  File "/usr/lib/python2.7/posixpath.py", line 98, in splitext
    return genericpath._splitext(p, sep, altsep, extsep)
  File "/usr/lib/python2.7/genericpath.py", line 99, in _splitext
    sepIndex = p.rfind(sep)
AttributeError: 'int' object has no attribute 'rfind'


texts: Creating WMS with QGIS Server using data from PostGIS

texts:
I am trying to publish a WMS via QGIS using QGIS Server. It works when the data for the QGIS project are shapefiles that I put in the same folder as the project on the QGIS Server (/usr/lib/cgi-bin/qgis-server/projectfolder).

Now I have data that are stored in a PostGIS database. When trying to load the WMS in OpenLayers it doesn´t load the requested WMS layer. I assume that the problem is that it fails to connect to the database since the password is not stored anywhere. I stored the connection to the database in QGIS Desktop but it looks like this information is not passed to the Server.

I am using:

QGIS-Version: 2.14.8-Essen
Apache-Version: 2.4.7
PostgreSQL-Version: 9.3.16
PostGIS: 2.1.2 r12389


texts: Connecting column from Excel table with attribute table in QGIS

texts: Is it possible to connect an Excel table or column from an Excel table with an attribute table in QGIS?


texts: List all grid squares crossed by multiple vectors

texts:
I have a vector table of road sections across the UK. I also have a vector table of 1km grid squares. For each road section, I want to report all the grid crossed. In QGIS Using the 'join attributes by location' function, it only reports the central grid square. However, some sections of road cross several sections and I want to report all squares crossed by the road section.

I am relatively newb to QGIS.


texts: How to use SLD to filter a WMS layer but without setting any styles?

texts:
I'm trying to filter (using OpenLayers) a WMS layer, served by MapServer.

I've got this working by building the appropriate OpenLayers.Filter rule(s) with an empty Polygon symbolizer (so no specific styles are set). I then construct an SLD from the rule(s) and pass that in the SLD_BODY parameter.

This has been working fine, however I've hit a snag. Originally, the map files for the various layers were setting STYLEITEM "AUTO", which was pretty ugly, but it worked.

I've now got rid of this and am manually setting the styles depending on the feature. e.g.

CLASSITEM "SomeField"
CLASSITEM "AnotherField"

CLASS
    EXPRESSION ("[SomeField]" eq "SomeValue" && "[AnotherField]" IN "val1,val2")

    STYLE
        COLOR 255 255 255
        OPACITY 25
    END #STYLE

    STYLE
        SYMBOL "hatch"
        COLOR 0 200 255
        ANGLE 45
        SIZE 2
        WIDTH 3
        GEOMTRANSFORM (buffer([shape],-1))
    END #STYLE

    STYLE
        OUTLINECOLOR 0 200 255
        WIDTH 1
    END #STYLE
END #CLASS

...Plus a good few more per layer.

Again, this all works fine. The layers look as I expect them too...until a filter is applied.

Then, the SLD_BODY gets generated with an appropriate filter and with no style data. However, my various styles defined (in the map files) then get ignored and every feature is just rendered in grey.

What am I doing wrong? It seems like the classes (or the expressions in the classes) are being ignored when using a filter.

The generated SLD document looks like:

<?xml version="1.0" encoding="UTF-8"?>
<sld:StyledLayerDescriptor xmlns:sld="http://www.opengis.net/sld" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
   <sld:NamedLayer>
      <sld:Name>MyLayerName</sld:Name>
      <sld:UserStyle>
         <sld:FeatureTypeStyle>
            <sld:Rule>
               <ogc:Filter>
                  <ogc:PropertyIsEqualTo matchCase="false">
                     <ogc:PropertyName>SomeField</ogc:PropertyName>
                     <ogc:Literal>SomeValue</ogc:Literal>
                  </ogc:PropertyIsEqualTo>
               </ogc:Filter>
               <sld:PolygonSymbolizer>
                  <sld:Fill />
                  <sld:Stroke />
               </sld:PolygonSymbolizer>
            </sld:Rule>
         </sld:FeatureTypeStyle>
      </sld:UserStyle>
   </sld:NamedLayer>
</sld:StyledLayerDescriptor>

And the appropriate part of the OpenLayers JS:

rule = new OpenLayers.Rule({
    filter:     layerFilter,
    symbolizer: {'Polygon': {
    }}
});

Side note:
If I manually change the <sld:Fill /> to <sld:Fill></sld:Fill>, it works. :|

Although I could just use JS to replace it, that seems like a really hacky, messy way of doing it and would be a real last resort.

What's the proper way of doing this?

Screenshots attached to show what I mean regarding styling:

Styles working (No filters)

Styles broken (Filter active for the main feature in the image)


texts:
when trying to activate GniPlannerFTTH plugin on QGIS, I always get an error calles

  "Couldn't load plugin GniPlannerFTTH due to an error when calling its classFactory() method"

and further

  ImportError: No module named openpyxl

I successfully installed scipy and openpyxl and there were no errors that have been returned when importing them in python. I've done anything that @SIGIS wrote in Couldn't load plugin GniPlannerFTTH.

In QGIS I also get the following traceback information and the Python path:

  Traceback (most recent call last):

  File "C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line
  333, in startPlugin plugins[packageName] = package.classFactory(iface)

  File
  "C:/Users/user/.qgis2/python/plugins\GniPlannerFTTH__init__.py", line
  31, in classFactory from gniplanner import GniPlannerFTTHPlugin

  File "C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line
  607, in _import mod = _builtin_import(name, globals, locals, fromlist,
  level)

  File
  "C:/Users/user/.qgis2/python/plugins\GniPlannerFTTH\gniplanner.py",
  line 34, in from gniplanner_provider import GniPlannerFTTHProvider

  File "C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line
  607, in _import mod = _builtin_import(name, globals, locals, fromlist,
  level)

  File
  "C:/Users/user/.qgis2/python/plugins\GniPlannerFTTHgniplanner_provider.py",
  line 32, in from planner_report import PlannerReport

  File "C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line
  607, in _import mod = _builtin_import(name, globals, locals, fromlist,
  level)

  File
  "C:/Users/user/.qgis2/python/plugins\GniPlannerFTTH\planner_report.py",
  line 26, in from openpyxl import Workbook

  File "C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line
  607, in _import mod = _builtin_import(name, globals, locals, fromlist,
  level)

  ImportError: No module named openpyxl

Python Path:

  C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python/plugins\processing

  C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python

  C:/Users/user/.qgis2/python

  C:/Users/user/.qgis2/python/plugins

  C:/PROGRA~2/QGIS2~1.18/apps/qgis/./python/plugins

  C:\PROGRA~2\QGIS2~1.18\bin\python27.zip

  C:\PROGRA~2\QGIS2~1.18\apps\Python27\DLLs

  C:\PROGRA~2\QGIS2~1.18\apps\Python27\lib

  C:\PROGRA~2\QGIS2~1.18\apps\Python27\lib\plat-win

  C:\PROGRA~2\QGIS2~1.18\apps\Python27\lib\lib-tk

  C:\PROGRA~2\QGIS2~1.18\bin

  C:\PROGRA~2\QGIS2~1.18\apps\Python27

  c:\progra~2\qgis2~1.18\apps\python27\lib\site-packages\jinja2-2.7.2-py2.7.egg

  c:\progra~2\qgis2~1.18\apps\python27\lib\site-packages\markupsafe-0.23-py2.7-win32.egg

  c:\progra~2\qgis2~1.18\apps\python27\lib\site-packages\shapely-1.2.18-py2.7-win32.egg

  C:\PROGRA~2\QGIS2~1.18\apps\Python27\lib\site-packages

  C:\PROGRA~2\QGIS2~1.18\apps\Python27\lib\site-packages\win32

  C:\PROGRA~2\QGIS2~1.18\apps\Python27\lib\site-packages\win32\lib

  C:\PROGRA~2\QGIS2~1.18\apps\Python27\lib\site-packages\Pythonwin

  C:\PROGRA~2\QGIS2~1.18\apps\Python27\lib\site-packages\wx-2.8-msw-unicode

  C:/Users/user/.qgis2//python

  C:\Users\user.qgis2\python\plugins\mmqgis/forms

  C:/Users/user/Desktop/FttX Planner

Im struggling to interprete these information. Might there be a problem with the directory?


texts: CMV - How can I zoom in further using ArcGIS online basemaps?

texts: I need to be able to zoom in further in my cmv web map but I can't figure out where this is set. Any suggestions on where to look? I am also wondering if it is possible for opacity in these basemaps.


texts: WMS - Multiple paramters in CQL_FILTER not working?

texts:
I'm doing a WMS request with 2 cql_filters. This is the resultant query:

  ?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=tmp_vista_expedientes&cql_filter=strSubstring(FecRegistro%2C0%2C4)%3D2018;strSubstring(codigo,0,3)=%27201%27&time=1573127644387&WIDTH=256&HEIGHT=256&CRS=EPSG%3A25830&STYLES=&BBOX=0
  ... 0

Basically I have 2 filters:

strSubstring(FecRegistro,0,4)=2018

And

strSubstring(codigo,0,3)='201'

Separating them with a semicolon ;. But only the first filter is applied.

  How can I make both filters be applied?


texts: Seeking Python API to assess accessibility to foodbank site?

texts:
The aim is to look at accessibility or lack thereof, in certain areas of the city, in order to inform future foodbank services.

I have been working in Python/Spyder, and using OpenRouteService API to calculate the routes.

However the project coordinator from the Foodbank, has recently suggested that it would be beneficial to add another layer of granularity to the project - using public transport data and timetables - as foodbank sessions run at distinct times through the week.

Can anyone recommend an appropriate API I could use in Python?

I have looked at https://developer.here.com/ which seems to be only available in other languages.

I am using Windows 10.


texts: Returning raw address from pandas column using Geopy?

texts:
I have a Pandas dataframe, and I'm calling a geopy function on each row of 'address_original' column to get as much detail as possible about every address

from geopy.geocoders import GoogleV3

geolocator = GoogleV3(api_key='', timeout=5,
                      )

df['full_address_geocoded'] = df['address_original'].progress_apply(geolocator.geocode())
print df['full_address_geocoded'][:5]

0    (Veale Rd, New Plymouth, New Zealand, (-39.095...
1    (Veale Rd, Delaware, USA, (39.8036422, -75.494...
2    (1068 Clearwater Valley Rd, Clearwater, BC V0E...
3    (1605 Pine St W, Stillwater, MN 55082, USA, (4...
4                                                 None

The issue is that I need more information than there is included when using geocode method. Geopy also has .raw function, which works perfectly on a single string and outputs the data about what region the address belong to, area, etc.

places = '1238 Davie St, Vancouver, BC'

geolocator = GoogleV3(api_key='AIzaSyBlNIvZTk-BpWDeX1FFXPbx6QwbNzZL80w')
location = geolocator.geocode(places, language='en')
print location.raw

{u'geometry': {u'location_type': u'ROOFTOP', u'bounds': {u'northeast': {u'lat': 50.6539919, u'lng': -120.3383232}, u'southwest': {u'lat': 50.6538198, u'lng': -120.3386968}}, u'viewport': {u'northeast': {u'lat': 50.6552548302915, u'lng': -120.3371610197085}, u'southwest': {u'lat': 50.6525568697085, u'lng': -120.3398589802915}}, u'location': {u'lat': 50.6539239, u'lng': -120.3385242}}, u'address_components': [{u'long_name': u'142', u'types': [u'street_number'], u'short_name': u'142'}, {u'long_name': u'Waddington Drive', u'types': [u'route'], u'short_name': u'Waddington Dr'}, {u'long_name': u'Upper Sahali', u'types': [u'neighborhood', u'political'], u'short_name': u'Upper Sahali'}, {u'long_name': u'Kamloops', u'types': [u'locality', u'political'], u'short_name': u'Kamloops'}, {u'long_name': u'Thompson-Nicola', u'types': [u'administrative_area_level_2', u'political'], u'short_name': u'Thompson-Nicola'}, {u'long_name': u'British Columbia', u'types': [u'administrative_area_level_1', u'political'], u'short_name': u'BC'}, {u'long_name': u'Canada', u'types': [u'country', u'political'], u'short_name': u'CA'}, {u'long_name': u'V2E 1N3', u'types': [u'postal_code'], u'short_name': u'V2E 1N3'}], u'place_id': u'ChIJaV6wtDgsflMR1J9ReM0lzLs', u'formatted_address': u'142 Waddington Dr, Kamloops, BC V2E 1N3, Canada', u'types': [u'premise']}

The raw method doesn't work on multiple values, however, I get the TypeError: 'dict' object is not callable.

What is a way around this?


texts: How to construct a line from points in pyshp?

texts: I am reading a shapefile for a pipeline and what I've noticed is that if I just draw a line from one point to the next (in the order I read them), there will be points where the line terminates, and the next point is really a new line - but I can't see any way to tell where those points are...  How does that work?


texts: Replace NA's with 0 for raster data using R?

texts: I am new in R. How can I replace NA's with 0 for my raster data which have spatial information?


texts: Error: Objects in this class cannot be edited outside an edit session

texts:
I'm simply trying to modify the Default Editing Widget ESRI sample to show my extent and service and am noticing the following error in Firebug when trying to add a point:

{"error":{"code":400,"message":"Unable to complete operation.","details":["Unable to perform applyEdits operation."]}}

I also get the following errors in my logs on ArcServer:

a database error occurred
objects in this class cannot be edited outside an edit session

Any ideas on why my data is locking like this? Or how I can fix it?


texts: Where can I get jsde_sdk.jar for ArcSDE 10?

texts:
I had an environment working with no problems, with Geonetwork harvesting metadata from ArcSDE 9.3. It was done trhu the following libraries, copied from the ArcSDE installation to the Geonerwork lib dir:

jsde_sdk.jar
jpe_sdk.jar
concurrent.jar
icu4_3_2.jar

After upgrading to ArcSDE10, the harvesting stopped working.

It is certainly because I'm using *.jars for version 9.3, while having a v10 ArcSDE, but the new installation doesn't have those jars.

Does anyone knows where to get them?


texts: Getting a square mile bounding box anchored by coordinates?

texts:
Let's say I have a lat/lng coordinate pair.  I want to determine the coordinates that would create a square mile box.  The existing coordintae pair would be the top left of the box, and the coordinate pair I'm trying to calculate would be the lower right.

What's the best way to handle this?  A solution using GeoDjango would be great, but PostGIS or even a mathematical formula would be good too.


texts: How to create a startup prompt that populates definition queries on multiple feature classes

texts:
I would like to create an ArcGIS start up prompt or custom button that asks for a user's initials and then populates multiple feature class & layer definition queries.

More specifically, I would like to standardize my departments layer organization for new employees and would like to customize the definition queries & symbology of several layers in each employee's MXD.


texts: Using customized Coordinate System in ArcGIS Desktop?

texts:
I don't know so much about coordinate systems...

In my office we use to deal with spatial data coming from archaeological sites. Each site has its own x-y-z coordinate system (GCS). Three simple ortogonal cartesian axis.
In the last years we have been managing this spatial data through GIS software (ArcGIS), without using specific coordinate system (just leave it as "undefined")

I'd like to know if there exisits any GCS designed to deal with such datasets using simple cartesian orthogonal axis, without grid distortions of the typical GCS. In addition, I'd like to know if this system is suitable for using it in an online mapping application.

By the way, we manage 2D (ArcMap) and 3D (ArcScene) environments and work with "mm" as length base unit.

If such a thing doesn't exists, maybe someone knows how to create it.


texts: Is there a way to avoid the re-calculation of the same centroid without creating a new column?

texts:
All tables contain polygon geometries.

SELECT *
FROM table1, table2, table3, table4
WHERE
    st_contains(table2.geom, st_centroid(table1.geom)) AND
    st_contains(table3.geom, st_centroid(table1.geom)) AND
    st_contains(table4.geom, st_centroid(table1.geom)) AND

Will the centroid be calculated more than once or is postgis/postgresql able to resolve it. If not, is there a possibility to improve the performance without creating a centroid column?


texts: Can ogr2ogr "reverse clip" (or "clip out" or "erase" or "difference") one shapefile from another?

texts:
Is there a means in ogr2ogr to "clip out" one shapefile from another?

I know ogr2ogr -clipsrc B.shp out.shp A.shp will remove any part of A.shp that is NOT in B.shp.

Is there an ogr2ogr command to remove any part of A.shp that intersects with B.shp?

I have seen this called "reverse clip", "clip out", "erase" (ArcMap), or "difference".

Note: this question is about ogr2ogr, the command line tool, NOT the OGR library, which requires programming.


texts:
So, I don't really know how to word what I want to ask, but I took a DEM file and did used the raster calculator to select a certain elevation range. After that I reclassified it to separate  it into different levels of elevation.

How do I get the new reclassified layer to have attribute data of those levels of elevation?


texts: Seeking System Requirements of GME for ArcGIS Desktop?

texts:
I've been trying to install GME 0.7.4 without success. This version is designed to run with ArcGIS 10.3 but I have installed ArcGIS 10.4.

Does anyone know if ArcGIS 10.4 is compatible (or can be made to be compatible) with GME 0.7.4?


texts: QGIS map file to be viewed on both mac and pc platforms

texts:
Im teaching a class where some students have QGIS installed on their macs and I'm using a pc.  Although the file structure is the same mac reference files differently than the pc and I get data link errors when I open a students map project.

Any quick conversion methods?  I can do it a link at a time but that seems a bit slow to me.


texts: Trouble saving graduated color style in Spatialite db in QGIS

texts:
When I attempt to save a style with graduated colors for a vector layer in a Spatialite db, I receive the following error:

  Error looking for style. The query was logged

I had success saving single symbol and categorized color styles following the workflow shown here. How can I save my graduated color class in the spatialite db?

QGIS: 2.12.1. Spatialite: 4.1.1


texts: How do I move the last two characters of a column within the attribute table to another column?

texts:
I have data in a column that looks like the below example:

El Monte CA

I need to move the last 2 characters to a separate column. Which function in the field calculator will deliver those results?


texts: Using raster calculator in ModelBuilder to batch reclassify hundreds of rasters?

texts:
I'm attempting to use iterate by raster to feed all the files through a simple map algebra expression (I'm looking at data that indicates probability of salamander habitat) so I want to reclassify the raster data and get output that separates 0 probability of habitat vs. greater than 0 probability.

After many hours of geoprocessing, I just realized that all the habitat extents are identical.

# -----------------------------------------------------------------------
# reclassifymodel.txt
# Created on: 2013-12-02 05:29:53.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: reclassifymodel <SalamanderData>
# Description:
# This model reclassifies all TIF rasters to binary probability
# -----------------------------------------------------------------------

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Load required toolboxes
arcpy.ImportToolbox("Model Functions")

# Script arguments
SalamanderData = arcpy.GetParameterAsText(0)
if SalamanderData == '#' or not SalamanderData:
    SalamanderData = "C:\\Users\\mkcar_000\\Desktop\\SalamanderData" # provide a default       value if unspecified

# Local variables:
Raster = SalamanderData
reclassified = Raster
polygon = reclassified
polyselect = polygon
extent__Name_feat_shp = polyselect
pathname = Raster
Name = SalamanderData

# Process: Iterate Rasters
arcpy.IterateRasters_mb(SalamanderData, "", "TIF", "RECURSIVE")

# Process: Parse Path
arcpy.ParsePath_mb(Raster, "PATH")

# Process: Raster Calculator
arcpy.gp.RasterCalculator_sa("outraster = Con(\"%Raster%\" > 0, 1, 0)", reclassified)

# Process: Raster to Polygon
arcpy.RasterToPolygon_conversion(reclassified, polygon, "SIMPLIFY", "Value")

# Process: Select
arcpy.Select_analysis(polygon, polyselect, "\"grid_code\" = 1")

# Process: Copy Features
arcpy.CopyFeatures_management(polyselect, extent__Name_feat_shp, "", "0", "0", "0")


texts: How do I generate a "granule" index shapefile for ImageMosaic in GeoServer?

texts:
I am experiencing a problem with generating a granule index shapefile for the ImageMosaic plugin in GeoServer; there are no clear instructions to follow.

Is there any way that I can generate it automatically using a raster dataset? I have a huge raster dataset to serve as WMS.


texts: Spatial join for multiple shapefiles (loop)

texts:
I have the following two datasets:

polygon shapefile containing administrative boundaries of a City
multiple point shapefiles

What I want to do is to spatially join the point shapefiles with the polygon based on the match_option intersect using arcpy. So to say, I want to count the number of Features per Point shapefile in the polygon shapefile. I do not Need any fields of the Point shapefils, but only the join_count field.

This is the python script I have so far.

output = 'F:\\test\\spatialJoin'
pointShpfiles = listofShapefiles.split(";")
spatialJoin = administrativeUnits
for pointShp in pointShpfiles:
    spatialJoin = arcpy.SpatialJoin_analysis(spatialJoin, pointShp, arcpy.Polygon, "JOIN_ONE_TO_ONE", "KEEP_ALL", "dummyField", "INTERSECT")
    arcpy.DeleteField_management(spatialJoin, 'dummyField')
    arcpy.DeleteField_management(spatialJoin, 'TARGET_FID')

arcpy.CopyFeatures_management(spatialJoin, output)

The script works partly, but the Problem is, that only the last of the Point shapefiles in the list is added, so I only have one join_count field instead of 3 (for three entered Point shapefiles).

What I am doing wrong?


texts: Problems merging OS Terrain50 tiles i.e. gaps at small scales after using "Mosaic to New Raster" Any ideas?

texts:
I downloaded the OS Terrain 50 dataset in ascii format, then unzipped all the folders and searched for all *.asc files.  I then copied and pasted these into a new folder, so that they could be easily located when geoprocessing.

I then used the "Mosaic to new raster tool" in ArcMap 10.2 and selected all the asciis as the input rasters, saved to a new geodatabase, defined projection (BNG), pixel type (32 bit float) and number of bands (1) to match the original files.  I left everything else as default options.

The resulting image has large areas missing in the south east of England (see below).  When you zoom in beyond 1:750,000 the gaps disappear.  Another strange thing is that when I removed the layer from the table of contents, then re-added it the whole dataset now only shows when zoomed in beyond 1:1,500,000 (and gaps still until 1:750,000).

I don't desperately need to show the data at national scale, however just curious as to what might be wrong and whether anyone can help.  Thanks, Sarah


texts:
I have a shapefile having 1000s of non overlapping polygons. I have got seamless layer of DEMs which cover the full extend of this polygon shapefile.

I want to find out the maximum altitude value in each of these polygons automatically. Is there any way to do in Qgis, ArcGIS, gdal or any other open source application.


texts:
I wrote some code for removal of polygon holes whose area value is under the hardcoded threshold.

source_g.getFeatures()[0].getGeometry().getPolygons().forEach(function(poly){

    if(!(poly.getLinearRings().length<=1))
    {
        poly.getLinearRings().forEach(function(ring){
            if(Math.abs(ring.getArea()) < 0.1){

               ring.setCoordinates([])
            }
        })
    };
})

The idea is to set coordinates of inner LinearRings to [], I tried to do that on the whole feature and it disappeared. But for LinearRings it's not doing the same.


texts: Adding unique IDs to mutiple polygons using loop or apply function in R

texts:
I want to add a unique ID (row number) for each polygon in a SpatialPolygonsDataFrame. I am new at working with spatial data (and programming). The goal is to merge all polygons belonging to a single species.

This is what I run:

setwd("~/type1")
species <- 'speciesX'
list_shp <- list.files(path=species, pattern="*.shp", full.names = TRUE,
                   recursive = TRUE, include.dirs = FALSE)
shp_objects <- lapply(list_shp, function(x) {readOGR(dsn=x, layer=ogrListLayers(x))})

# See the row names (IDs)
lapply(shp_objects, function(x) row.names(x@data))

# Manually set unique row names (IDs)
shp_objects[[1]] <- spChFIDs(shp_objects[[1]], as.character(10)) # e.g. 10

But, since the SPDF for each species can have several polygons, I tried making a for loop, implementing a simple counter, adding a consecutive number to each as row names (IDs) using the spChFIDs function.

idx<-0
for(i in shp_objects){
 idx<-idx+1
 shp_objects <- spChFIDs(i, as.character(idx))
}

When running this, and then lapply(shp_objects, function(x) row.names(x@data)) to see if the row names have been updated accordingly, I get the error:

Error in as.list.default(X) :
  no method for coercing this S4 class to a vector

How could I make this work?


texts:
I've got a map laid out how I want it in QGIS 2.2, using the "Atlas Generation" feature. When I preview the atlas, all 7 pages look proper and seem to work fine. However, I cannot export them to any format as it always results in blank pages (the proper number though, so it appears to be 'iterating' properly).

If it matters, it's set to print as raster as I have some transparency.  Thoughts on what is wrong?


texts: QGIS 3 heatmap, on screen color ramp with data scale

texts:
When making heatmaps in qgis 3.0.1 i would like to place the color ramp with percentage and actual numbers for the bin borders on the map.

Each marker under the color ramp represent a percentage of the max value. So at a maxvalue of 4000, the marker at 50% would represent 2000.

Is it posible to access the actual number for each marker/bin or at least the maximum value for the heatmap/color ramp? One could then do the rest manualy.


texts:
I have download today the Google Earth engine plugin to QGIS using  this:
https://gee-community.github.io/qgis-earthengine-plugin/

I have tried to upload the image as shown in the example :

I have coped the code but for some reason I get erro and can't add the layer:

import ee
from ee_plugin import Map
image=ee.ImageCollection("COPERNICUS/S2").filterBounds(Map.getCenter()).first()
Map.addLayer(image, {'min':500, 'max':5000, 'gamma': 1.4, 'bands': ['B4','B3','B2']}, 'image')

then I get the error:

  raise ee_exception.EEException(json_content['error']['message'])
  ee.ee_exception.EEException: Image.visualize: Parameter 'image' is
  required.

I believe I have define the visualization and that there is an image, so why do I get this error?

My endgoal is to be able to upload images to QGIS using GEE plugin


texts: CSV list as input parameter in QGIS 3.0 Processing Modeler

texts:
I have CSV list with XY coordinates (company and its competitors) which I want to input into simple QGIS Processing model. From this list the model creates and styles two spatial layers.

Everything works well with one exception. I have the input CSV list set up to fixed directory path (see screenshot no. 1). I would like if the model allowed me to select the path to the CSV list when I run the model (see screenshot no. 2).

Does anyone know how to set up a CSV list as an input parameter to QGIS processing model?


texts: Filtering GPS data for precise positioning?

texts:
I am using GPS receiver named G STAR-IV Globalsat BU-353-S4 but it gives erroneous position.

How can I get better position in Real Time using any filtering techniques?

I saw some filtering algorithms but most of the algorithms talk about base station and I don't know about base station. Please see the image. The green is the true path and red points are GPS position data.


texts: How to find the furthest coordinates on a polygon in a given direction

texts:
My task is to generate a box that maps over the existing polygon, which I’ll call the level. What I’m doing is mapping buildings with geojson, which does so by using a set of lat/long coordinates as points to define the polygon that matches the walls.
To create the box, I currently have a script in place to find the furthest coordinates (on the level) in North/South/east/West directions. It then generates a line with those two coordinates and finds the intersection between the line in each direction. Essentially boxing in the entire level no matter what funky shape the building is in.
However, this only works when the building directly faces a cardinal direction. My question is, when the building is instead facing X degrees in a certain direction, how can I perform the same calculation that will generate the box for me?


texts: Fastest methods for modifying attribute tables with Python?

texts:
A while ago, I wrote a quick Python function for converting an attribute table to a python dictionary, where the key is taken from a user-specified unique ID field (typically the OID field). Additionally, by default all fields are copied to the dictionary, but I've included a parameter allowing for just a subset to be specified.

def make_attribute_dict(fc, key_field, attr_list=['*']):
    dict = {}
    fc_field_objects = arcpy.ListFields(fc)
    fc_fields = [field.name for field in fc_field_objects if field.type != 'Geometry']
    if attr_list == ['*']:
        valid_fields = fc_fields
    else:
        valid_fields = [field for field in attr_list if field in fc_fields]
    if key_field not in valid_fields:
        cursor_fields = valid_fields + [key_field]
    else:
        cursor_fields = valid_fields
    with arcpy.da.SearchCursor(fc, cursor_fields) as cursor:
        for row in cursor:
            key = row[cursor_fields.index(key_field)]
            subdict = {}
            for field in valid_fields:
                subdict[field] = row[cursor_fields.index(field)]
            dict[key] = subdict
            del subdict
    return dict

This works great for relatively small datasets, but I just ran it on a table containing about 750,000 rows and 15 fields -- around 100MB in a file geodatabase. On these, the function runs much slower than I would have expected: around 5-6 minutes (and this is after copying the table to the in_memory workspace). I'd really like to find a way to speed up the conversion to dictionary, or get some insight on a better strategy for manipulating large amounts of attribute data using Python.

UpdateCursors won't work well for me, because when one row changes, it has the potential to trigger changes in several others. Looping through and processing them one at a time is too cumbersome for what I need.


texts: QGIS SimpleSymbolRenderer label for Symbol in Legend

texts: Is there an ability to set a label for the symbol (in QgsSingleSymbolRendererV2) which is showing in the legend?


texts: Setting coordinate systems in QGIS: "transform error caught"

texts:
I'm learning how to use QGIS (v 1.8.0 "Lisboa") on Mac OS 10.6.8. My goal is to create a simple county map of Illinois using the base data shapefiles available from the Illinois Natural Resources Geospatial Data Clearinghouse.

Following workflow suggestions in another thread, I first checked "Enable 'on the fly' CRS transformation" in Project Properties.

Then I added the shapefile:

I then changed the CRS in Project Properties to NAD83/Illinois (West) [EPSG code 26972] to get the desired look:

So far so good. But after saving the QGIS file, closing it, then opening it again, I get a CRS error:

Transform error caught: forward transform of (0.000000, 0.000000)
PROJ.4: +proj=longlat +datum=WGS84 +no_defs +to +proj=tmerc +lat_0=36.66666666666666
+lon_0=-90.16666666666667 +k=0.999941177 +x_0=700000 +y_0=0 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Error: latitude or longitude exceeded limits

How do I prevent this type of error?


texts:
I've been given a number of files including some annotated shapefiles - Basically kilometre post markers and cable identification markers, to go with a polyline shapefile.  I was hoping that we could just drag 'n' drop the files over the images and lines in QGIS map window and they would show, but we get an error message (unknown file type) or the annotations just show up as blobs.

The files are .DBF files, .PRJ files, .SBN files, .SBX file with an autoCAD LT Shape source and AutoCAD LT Compiled Shape, together with an .XML file.

how can I overlay the annotations and get them to show text rather than the blob?

I can click on the properties in the layer browser, and see the text that is supposed to display, but I can't work out how to get it to show.


texts: Loop through 2 folders and perform raster calculation

texts:
I am currently learning arcgis and python.

The problem includes looping through 2 folders (with rasters of same filename) and performing raster calculation on each raster pair. Below is my code:

import arcpy, os, glob
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")

ws1 = "E:\clip_qx"
ws2 = "E:\clip_qy"

outws = "E:\result"

rasterlist = arcpy.ListRasters()
for r in ws1:
    r1 = arcpy.sa.Raster(r)
    r2 = arcpy.sa.Raster(os.path.join(ws2, basename))
    result = sqrt ( r1^2 + r2^2 )
    outname = basename
    result.save(os.path.join(outws, outname))

However, I am getting an error. It seems that the rasters were not read.

RuntimeError: ERROR 000732: Input Raster: Dataset E does not exist or is not supported.


texts: Finding value of peer to peer pixel in another raster layer with same pixel size, coordinate and location from attribute table?

texts:
I have two raster layers with same pixel size, location, coordinate and time. I need to know what the value of defined pixel in one layer is in same location (pixel) in another raster layer.
The manual way in ArcGIS is i (identify) ,however there are 144 layers! I need an easier way

I work with ArcGIS for Desktop.

The images are LST and MODIS NDVI. The image of NDVI was Re-sampled.The pixel size is 5000m. I need to know the value of the right image.

[![enter image description here][1]][1]


texts:
I tried many times to create new layer from query result, but i got this error message, when loading it.

" PostgreSQL layer has no primary key."

    SELECT * FROM pgr_drivingDistance('select gid::int4 as id, source, target,
length_m as cost from public.ospr_ways', 24710, 1000
) AS routed1
INNER JOIN ospr_ways AS e
ON routed1.edge = e.gid::int4;

update: i have installed QGIS 3.10.5 LTR 64x, PostgreSQL 12 64x, PostGIS 3.0 for PostgreSQL 64X, Windows 10 64x
i already created postgis, pgrouting extension from pgAdmin for target database, and the query successfully executed and displayed it from pgAdmin, but i can't load my query as a new layer from db manager of QGIS.


texts: QGIS Raster calculation with nodata values?

texts:
I am using QGIS 2.18.13.

I have a few shapefiles. From inside QGIS, using python's processing.runalg function and supplying "grass7:v.to.rast.attribute", I converted those to raster files which have value of 1 everywhere shapes were present and "no data" (255) beyond them.

Those rasters share the same spatial location and I would like to "add them" so if in the same location there are 3 values 1 in final raster there will be value 3. So, I will be happy to have nodata converted to 0s for that purpose.

I tried following things where I didn't get results, I found some limitations or I got no results at all.

Raster -> raster calculator - I don't have any control over
resolution. It was produced uncompressed tiff (400MB! from my
10000x10000 raster), 32 bit precision and nodata values changed
to -3.40282e+38. Where it was nodata before I get nodata so,
resulting raster is some kind of intersection rather than sum.
Toolbox -> GRASS GIS 7 -> Raster -> r.null. I get this error: "The following layers were not correctly generated. NullRaster"
toolbox -> GRASS GIS 7 -> Raster -> r.mapcalc. I don't know how to write an expression. I always get either parse error or:

r.mapcalc expression="abc = 1" --overwrite
r was unexpected at this time.

there are usually 4 lines in the log like this:

ERROR: Sorry is not a valid option

toolbox -> SAGA -> Raster calculus -> Raster calculator - that one seems not to throw errors, however, I am not sure how to handle NaN's?

I would like to add this to python script so I can automate the process if my base datasets change.


texts:
Does anyone know of an API (preferably free) that can provide real time world statistics such as Crude Birth Rate, Crude Death Rate, world/country population, or emails sent today as on the website http://www.worldometers.info/?

I'm looking to connect the data to an iPhone app, and I can work with any sort of return type.