texts: Are there any continental-scale flow accumulation products besides Hydrosheds?

texts:
I am trying to compute branching angles of drainage networks for roughly the western half of the U.S. Because many of these streams are non-perennial (don't have water all year), some of the standard data products won't work (Openstreetmaps, National Hydrology Dataset, Global Surface Water rasters). It seems the best approach is to use a flow accumulation raster.

I am aware of the Hydrosheds 15- and 30-arc second global flow accumulation grids, but they are too coarse. The higher the resolution, the better, but it looks like 10m DEM from the National Elevation Dataset is the best I can do. I am wondering if I can skip the steps of procuring and processing the data as I only need the flow accumulation. Deriving flow accumulation grids from DEMs seems like such a common task that I thought there would be more continental-scale flow accumulation grids available, but I can only find the Hydrosheds. Does anyone know of larger-scale products?


texts: Merge intersecting geometries from 2 tables into third table

texts:
I am attempting to take data from 2 tables, find intersecting geometries, and insert rows into the new table representing the intersections.  I have done this in the past with CartoDB using SQL but now I'm not able to do so.  I'd also be open to using the Merge Dataset functionality via the UI, but that doesn't appear to be designed for what I need.

Specifically, when I run the following sql via the UI, I get an error stating "SELECT ... INTO is not allowed here".

SELECT clipped.city_name, the_geom
INTO newtable
FROM
    (SELECT table1.city_name, table2.namelsad, (ST_Dump(ST_Intersection(table2.the_geom, table1.the_geom))).geom As the_geom
    FROM table2
    INNER JOIN table1
    ON ST_Intersects(us_county.the_geom, us_or_city.the_geom))  As clipped

Any suggestions on how to make this work?


texts: Why the georeferencing in UTM does not work with the QGIS georeferencer?

texts:
First time here, sorry if this looks simple.

I use QGIS on Windows 7. I georeferenced a toposheet with geographic references to UTM. When I superimpose a UTM raster, the toposheet is near the origin and the rater far away but properly imported.

I then georeferenced the toposheet with its system and warped it to UTM and it worked, the two rasters are properly superimposed. Can someone explain why the georeferencing in UTM does not work the first time with the georeferencer?


texts: Listening to ICommandBar or ICommand OnClick event?

texts:
I have an ArcMap extension and I'm listening to the IEditEvents.

When the OnChangeFeature, OnCreateFeature or OnDeleteFeature event is called, I want to know if it came from the "merge(Editor_Merge)" button.

Is there a way I can listen to the onClick event of the merge button?


texts: Adding timestamp to video export in Google Earth Engine

texts:
I am exporting an image collection to video using:

Export.video.toDrive({
  collection: mov_EVI,
  description: 'Modis EVI-Classified Water',
  dimensions: 720,
  framesPerSecond: 10,
  region:ROI
});

How can I add a timestamp to the export frames of the video using system:time_start of each image?


texts: Changing polygon boundary thickness/color

texts: How can I make the polygon boundaries of admin 1 thicker than thous admin 2? I. e. I want to make the province boundaries thicker than its district boundaries.


texts: Use MapBox Satellite basemap in ArcGIS Online web map

texts:
I want to add the MapBox Satellite basemap layer as an option in my ArcGIS Online web app.  There is a widget that allows users to toggle between basemaps and I would like MapBox Satellite to be an option there.  Is this possible?  I know the URL JOSM is using to display it is this:

tms[19]:http://{switch:a,b,c}.tiles.mapbox.com/v4/openstreetmap.map-inh7ifmo/{zoom}/{x}/{y}.png?access_token=pk.eyJ1Ijoib3BlbnN0cmVldG1hcCIsImEiOiJncjlmd0t3In0.DmZsIeOW-3x-C5eX-wAqTw

Is there anyway to tweak this so I can use it within the basemap widget?
Here is the map I want to use it in.


texts: Overlapping polygons and features and how to dissolve

texts:
I am new to GIS area and also new to any software that can modify the GIS.

I am currently working with South Australia map I have download from Australian census site. I am using QGIS to load the shapefile and also do any modification.

Here is a link to download the geo datasource from Australian Government Site

http://data.gov.au/dataset/bcfcfc9a-7c8d-479a-9bdf-b95ca66ad29a/resource/5ce6c946-787b-4c4e-ba6f-c2ec18dfa415/download/salocalitypolygon.zip

My problem is when I have downloaded this south australia (SA - State of Australia), near the boundary of the state, I see there two polygons on top of each other fully overlapping. I reduced the fill opacity to 50% so where you see dark yellow colour those are overlapping polygons. It looks like some rectangle polygons on top of real suburb's outline. I dont know why it is like that ? And how can I get read of this.

Also I want to know how can merge (or dissolve) features (polygons) in same layer. Everytime I dissolve it create on new layer.


texts: Settings in ArcMap for exporting maps to be used in PowerPoint?

texts:
I was wondering if there are any good rules which can be adopted when making maps that are only going to be used digitally in PowerPoint, and what the best format to export is (tiffs, pdfs, jpegs)?

Does anyone have tips? I guess this will depends on detail, size, and aspect of what you are trying to show, but I was wondering if there are general rules or guidelines?


texts: How to use comparison and logical operators with filters for a SQL search in Google Earth Engine in a single line

texts:
Recently, I asked how to query by using an expression with comparison and logical operators (u'"area" > 2e8 AND "area" < 3e8') in PyQGIS 3. Searching for an equivalent expression in Google Earth Engine (GEE), I found out that it can be done with filters 'gt' and 'lt'. However, for obtaining same result I had to split expression in two parts:

// Filter to get "area" > 2e8 AND "area" < 3e8
var selected1 = polygon8.filter(ee.Filter.gt('area', 2e8));
var selected2 = selected1.filter(ee.Filter.lt('area', 3e8));

as it can be observed at following code:

// Load FeatureCollection

var polygon8 = ee.FeatureCollection("users/joseguerreroa/vector_files/polygon8");

// Center the map on the FeatureCollection.
Map.centerObject(polygon8, 9);

// Filter to get "area" > 2e8 AND "area" < 3e8
var selected1 = polygon8.filter(ee.Filter.gt('area', 2e8));
var selected2 = selected1.filter(ee.Filter.lt('area', 3e8));

// Display the FeatureCollection and print all queried elements.
Map.addLayer(selected2, {color: 'blue'}, 'polygon8');

When I use this line in my code:

var selected2 = polygon8.filter( ee.Filter.gt('area', 2e8) && ee.Filter.lt('area', 3e8));

code works but result is different as expected. So, my question is how to use comparison and logical operators with filters for a SQL search in Google Earth Engine in only a single line.

Editing note:

With 'filter', instead 'filterMetadata', it also works:

// Load FeatureCollection

var polygon8 = ee.FeatureCollection("users/joseguerreroa/vector_files/polygon8");

// Center the map on the FeatureCollection.
Map.centerObject(polygon8, 9);

// Filter to get "area" > 2e8 AND "area" < 3e8
var selected = polygon8.filter(ee.Filter.gt('area', 2e8)).filter(ee.Filter.lt('area', 3e8));

// Display the FeatureCollection and print all queried elements.
Map.addLayer(selected, {color: 'blue'}, 'polygon8');


texts: Extracting metadata from shapefile in FME?

texts:
Does anyone know how to extract metadata information from a shapefile in FME?

To be exact, i want to extract the schema (data structure) of the data such as the field names, etc. from a shapefile and export as a xml or any format.

Example, if i have a shapefile named 'Road', i want to extract its field names like road name, road length etc and export it in xml preferably rather than viewing the field names of the shapefile in a datainspector.


texts: postprocessing GPS coordinates with RINEX data to improve precision?

texts:
Im doing some mapping for OSM in my region and I am curious if the precision of tracks from a GPS Receiver can be improved if they are post processed with RINEX data.

The tracking is done with a Garmin GPSMAP 60CSX or alternatively my Android Phone and OSMTracker for Android.

The RINEX data is available from http://igs.bkg.bund.de/filestatus/dailyrinex

cheerio
   steve


texts: Linear referencing in Geographic Coordinate System (GCS) using ArcGIS for Desktop?

texts:
I want to make a new line segment in event table with linear referencing ArcGIS 10.2 for Desktop.

I'm using two tools:

Create Route with input line feature in Geographic coordinates system and tool
Make Route Event Layer with input table event with decimal degree unit.

But when I trying running them, the result is that only half of the events are in the segmentation.

Why are all lines not segmented based on table events?

in segment with low risk i give blue color and high risk in red color, but only low risk display in segment


texts: How can I select feature, except null values?

texts:
I want to select feature except NULL values.

In the past, when I using ArcGIS, I used 'select by attributes', and type like 'NOT "field name"=NULL'

But now, in QGIS, it doesn't work.
Just like ArcGIS, i used 'select by attributes', and type like 'NOT "field name"=NULL'
But it select nothing.

How can I select feature except NULL values, using 'select by attributes' in QGIS?


texts: Superposition of two cached maps with Android SDK

texts: I'm developping an Android Application with ArcGIS Runtime For Android, i wish to know if it's possible to superpose 2 cached maps from 2 diffrent WMS and store them in mobil.


texts:
In a couple of my scripts where I need to write a geolocalised raster using the GDAL library in Python, I need to run the script twice to obtain an output. The first run generates the file, but it is empty and comes up with a size of 0 octets on my linux system. It's only when I run the script a second time that the raster gets filled with the information. Why?

I am using anaconda with Python 3.4 and GDAL 2.0.0 and libgdal 2.1.1 (I had many issues finding the right versions for compatibility).

Below a working example that causes the issue for me:

import numpy as np
from osgeo import gdal, osr

# Create a random 100x100 raster
test_dem = np.random.rand(100,100)

# Get size of raster
DEM_size = np.shape(test_dem)

# Output file path
dstfile =  '/output/file.tiff'

# Make geotransform
xmin,ymax = [296058.21,4990799.17]
nrows,ncols = DEM_size
xres = 1
yres = 1
geotransform = (xmin,xres,0,ymax,0, -yres)

# Write output
driver = gdal.GetDriverByName('Gtiff')
dataset = driver.Create(dstfile, ncols, nrows, 1, gdal.GDT_Float32)
dataset.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(32632)
dataset.SetProjection(srs.ExportToWkt())
dataset.GetRasterBand(1).WriteArray(test_dem)


texts:
I have 2 tables, the table "fruit" contains a list of fruits (see figure 01).

The "municipality" table shows in the column 'execucao01' which fruits exist in each municipality (see figure 02).

I would like to create a virtual field in the "fruit" table to obtain the total existence of each type of fruit in the "municipality" table.
I tried to use the aggregation function.
count (aggregate (layer: = 'municipality', aggregate: = 'concatenate', expression: = "execucao01"), filter: = (array_contains (string_to_array (substr (aggregate (layer: = 'municipality', aggregate: = 'concatenate ', expression: = "execucao01"), 2, -1)),' "FRUIT" ')) = true)


texts: How can I come up with the spatial variation of ndvi using arcgis 9.3.

texts: I would like to come up with the spatial variation of NDVI of a certain area. I have the average NDVI from 1998 to 2014. I'm using Arcgis 9.3


texts: How do you create a "donut hole" / cut-out-effect using folium?

texts:
I've tried adjusting parameters and changing the order in which these circles are added. I can't get the interior circle, or "hole", to remain transparent. The outer circle always predominates.

m = folium.Map(tiles="cartodbpositron", location=xy, zoom_start=14, attr="CartoDB")

folium.Circle(xy, radius=1800, fill=True, fill_opacity=.5, fill_color='black').add_to(m)
folium.Circle(xy, radius=1000, fill=False, fill_opacity=0, fill_color=None).add_to(m)

m

In effect, I'm trying to do what Shapely does:

donut = Point(0, 0).buffer(2.0).difference(Point(0, 0).buffer(1.0))

But I want to do it with Folium because projecting a circle made with Shapely yields an object that's rough, even elliptical.


texts:
I am trying to add a feature layer to my map document as the output from my python script.  The goal is to limit the number of fields in the attribute table to a pre-defined list.  I am able to add the feature layer to my map document, but all of the fields display. However, when I use arcpy.Describe() on the layer my script added to the map, the Describe results show the desired field subset.  I can't figure out what is going on.

def makeFeatureLayer(outputFC, fieldTable, input_FC, input_ID, mapLayers):
    outputFeatureLayer = "{0}".format(os.path.basename(outputFC))
    subField = "{0}".format(arcpy.Describe(outputFC).subtypeFieldName)
    displayFieldList = [input_ID, subField]
    with arcpy.da.SearchCursor(fieldTable, ["layer_name", "Attribute"]) as cursor:
        for row in cursor:
            if "{0}".format(row[0]) == "{0}".format(os.path.basename(input_FC)):
                displayFieldList.append("{0}".format(row[1]))

    outputFieldList = arcpy.ListFields(outputFC)
    fieldInfo = "{0}".format(input_ID) + " " + "{0}".format(input_ID) + " " +  "VISIBLE" + " " + "NONE;" + "{0}".format(subField) + " " + "{0}".format(subField) + " " +  "VISIBLE" + " " + "NONE;"
    for f in outputFieldList:
        if f.name in displayFieldList:
            fieldInfo = fieldInfo + f.name + " " + f.name + " " + "VISIBLE" + " " + "NONE;"
        else:
            fieldInfo = fieldInfo + f.name + " " + f.name + " " + "HIDDEN" + " " + "NONE;"

    arcpy.management.MakeFeatureLayer(outputFC, outputFeatureLayer, "", "", fieldInfo)

    mapLayers.append(outputFeatureLayer)

    del subField, displayFieldList, outputFieldList, fieldInfo

def updateMap(mapLayers):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    dataFrame = arcpy.mapping.ListDataFrames(mxd)[0]
    for layer in mapLayers:
        arcpy.mapping.AddLayer(dataFrame, arcpy.mapping.Layer(layer), "TOP")
        arcpy.RefreshActiveView()
        arcpy.RefreshTOC()
    del mxd, dataFrame

I tried leaving out arcpy.mapping.Layer in case this was undoing the fieldInfo settings/field selection, and just trying to add the feature layer to the map document directly, but that generates an assertion error.  Also, if I run the code it actually results in two new layers to my map document: a visible layer that shows all the fields, and the original output of the MakeFeatureLayer tool, but it is not visible and cannot be added to the map display unless I export it to a .lyr file first.

Any ideas?

I was able to get MakeFeatureLayer to work from the Python window just fine.  I tried to use ListLayers() and use that to get at the feature layer after creating it, but that didn't work either.

I can sort of work around this by saving the results of MakeFeatureLayer to a .lyr file, but this is not ideal as I do not want to write anything extra to disk just for the field display. When I add the .lyr file to the map document, though, it works the way I intend and only shows the subset of fields.


texts: Exporting shapefile into PostGIS via Python script from QGIS

texts:
I want to export a shapefile into PostGIS via a Python script from QGIS. I’m using QGIS 3.10, PostgreSQL 10 and PostGIS 2.5 on Windows 10. I tried go use qgis:importintopostgis:
import processing
processing.run("qgis:importintopostgis",{
                'INPUT': 'unioncrs',
                'DATABASE': "HistoPG",
                'SCHEMA': "public",
                'TABLENAME': "testunion",
                'PRIMARY_KEY': "Nummer"
                })

That produces an error message:

“.... There were errors executing the algorithm.”

Interestingly it works with the DB Manager, which seems to be a GUI of qgis:importintopostgis.
The differences which I can see are the Source and the Target SRIDs. I don’t know how to enter these parameters in qgis:importintopostgis. I’m also not sure if I wouldn’t need to enter username and pasword for the database, but again I don’t see where to enter these parameters in qgis:importintopostgis.
I also tried `gdal:importvectorintopostgisdatabaseavailableconnections:
processing.run("gdal:importvectorintopostgisdatabaseavailableconnections",{
                'DATABASE': "HistoPG",
                'INPUT': 'C:/QGIS/Vers1/Results/unioncrs.shp',
                'GTYPE': 0,
                'SCHEMA': "public",
                'TABLE': "testunion"
                })

That runs without error message but also doesn’t import a table into PostGIS.
Is there a method which works with Python from QGIS?
I run:
PARAMETERS = {'DATABASE':'HistoPG','INPUT':'C:\\QGIS\\Vers1\\Results\\unioncrs.shp',
'SHAPE_ENCODING':'utf-8','GTYPE':0,'A_SRS':QgsCoordinateReferenceSystem('EPSG:31467'),'T_SRS':None,'S_SRS':None,
    'SCHEMA':'public','TABLE':'testunion','PK':'Nummer','PRIMARY_KEY':'','GEOCOLUMN':'geom','DIM':0,'SIMPLIFY':'',
   'SEGMENTIZE':'','SPAT':None,'CLIP':False,'WHERE':'','GT':'','OVERWRITE':True,'APPEND':False,
    'ADDFIELDS':False,'LAUNDER':False,'INDEX':False,'SKIPFAILURES':False,'PROMOTETOMULTI':True,
    'PRECISION':True,'OPTIONS':''}

processing.run("gdal:importvectorintopostgisdatabaseavailableconnections", PARAMETERS)

I then get the following error message:

Traceback (most recent call last):   File
"C:\PROGRA~2\QGIS3~1.10\apps\Python37\lib\code.py", line 90, in
runcode
exec(code, self.locals)   File "", line 1, in    File "", line 98, in    File
"C:/PROGRA~2/QGIS3~1.10/apps/qgis-ltr/./python/plugins\processing\tools\general.py",
line 106, in run
return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)   File
"C:/PROGRA~2/QGIS3~1.10/apps/qgis-ltr/./python/plugins\processing\core\Processing.py",
line 181, in runAlgorithm
raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.

It creates a table in my PostGIS database. But the table is empty and has only the fields ‚nummer‘ and ‚geom‘.
Do I have to define the fields somewhere in the parameters or do I have to append to a table which has the fields already?
I found the History, that’s good to know. But it doesn’t log when I run DB Manager. And I didn’t find any other plugin to import a vector layer into PostGIS.
Which one did you refer to?
I found the function now in the Processing toolbox. I ran it successfully and used the text from the history in Python. That produced after a while an error message „...There were errors executing the algorithm.“ But it had exported the table into PostgreSQL. The number of features is correct and the table looks alright so far. So I’m not quite sure what went wrong actually.


texts: sometimes a white, square dialogue appears in the middle of the screen: Map Data ... Google maps

texts:
I have this code and some of the time,(not always!!), i have a dialogue that get

displayed in the middle of the map:

Google Maps

Map data (c)2012 Google

   var map;
   var vectors;
   var epsg900913 = new OpenLayers.Projection('EPSG:900913');
   var epsg4326 = new OpenLayers.Projection('EPSG:4326');

   function initialize() {
       var control, renderer=OpenLayers.Util.getParameters(window.location.href).renderer;
       renderer= (renderer) ? [renderer] : OpenLayers.Layer.Vector.prototype.renderers;
       // Create the map object
       map = new OpenLayers.Map('map');
        //Create a Google layer
       var gmap = new OpenLayers.Layer.Google(
           "Google Streets", // the default
             {numZoomLevels: 20}
       );
       var wms = new OpenLayers.Layer.WMS( "OpenLayers WMS",
                "http://vmap0.tiles.osgeo.org/wms/vmap0?", {layers: 'basic'});

       vectors= new OpenLayers.Layer.Vector(
          "Vector Layer",
           {
               renderers:renderer
           }
       );

       map.addLayers([gmap,vectors]);

       map.addControl(new OpenLayers.Control.LayerSwitcher());
       map.addControl(new OpenLayers.Control.MousePosition());

       control=new OpenLayers.Control.DrawFeature(vectors,OpenLayers.Handler.Path);
       map.addControl(control);
       control.activate();
       // Zoom to Vancouver, BC
       map.setCenter(new OpenLayers.LonLat(-123.12, 49.28).transform(epsg4326, epsg900913), 13);
  }

Does anyone know where is this coming from?

Thanks in advance.


texts: Maximum number of inputs that ArcGIS Merge tool will accept?

texts: I've done an internet search trying to find the theoretical maximum number of input datasets that the Merge tool will accept.  Nothing has turned up.  I'm currently using 10.5.1


texts: QGIS field calculator returning results in different units?

texts:
I am trying to use QGIS to get the perimeter, area, bounds width, and bounds height of US congressional districts, using the US Census TIGER shapefile. I loaded the shapefile into QGIS, set my CRS to NZGD49 with on-the-fly transformation enabled, and then started doing calculations using field calculator.

Perimeter and area seemed to work fine; I got pretty large numbers as you'd expect for units in meters and square meters. But when I tried to get bounds width and bounds height using bounds_width($geometry) and bounds_height($geometry), I got very small numbers, mostly less than 5 - it looks as if the units are in degrees for bounds but in meters for perimeter and area. I tried different projections, but with similar results for all of them. How can I get QGIS to report bounds width and bounds height in meters instead of degrees?


texts: elevation difference along the border of two categories GRASS

texts:
Consider the following raster with two categories and a DEM (numeric raster), I'd like to compute the elevation difference along the border (closest pixel) of those two categories, it doesn't matter on which side of the border the result is stored. I've checked r.neighbors but it's not category-aware. It's of little matter what it does in corners, where there are two closest pixels.

I'd like to keep it in GRASS (maybe python) because I'm making lots of spatial stats through rgrass7.


texts: Comparing field values for soil polygons using ArcPy?

texts:
Occasionally NRCS will update their soils layer for each county.  When they do, we get an email from our clients because they have been notified of the change.  The problem is, sometimes it's just a metadata change.  And we don't know until we run the full analysis and compare the soil types (MUSYM).

I'm trying to write a script that will compare the soils layer on file with the supposed "new" soils layer. If the unique values in the MUSYM field are the same for both features, I would like it to print "no new soils".  If not, print "new soils added".  The data is stored in a File Geodatabase.  I'm using version 10.5 of Arcmap.  I've looked through some of the other posts but I'm too green yet to figure this out.

import arcpy, os
arcpy.env.workspace = "H:\\SoilsUpdate\\SoilsUpdate.gdb"

#INPUTS:
fc1 = "R:\\State Soils\\gSSURGO_KS.gdb\\MUPOLYGON"

fc2 = "C:\\County\\Assessor.gdb\\NRCS\\Soils"

#set up cursors
cursor1 = arcpy.da.SearchCursor(fc1, ["MUSYM"])

cursor2 = arcpy.da.SearchCursor(fc2, ["MUSYM"])

#make a dictionary and store values
soil_lst = {}

for row in cursor1:
  soil_lst[row[0]] = row[1]

#loop through other table
for row in cursor2:

  if row[1] = soil_lst[row[0]]:

   soil_lst.append(row[0])

    print "no new soils".format(row[0])


texts: LIKE operator is not case-sensitive in Filter query for shapefiles (QGIS 3.4)

texts:
When I try to filter local shapefile with Query Builder (right-click on layer > Filter) using LIKE operator (which should be case-sensitive) it does not work properly. For example
query "place_name" LIKE '%Wola' returns both Stara Wola and Starowola (identically to ILIKE).
I don't know why, but it works correctly when the same filter query is performed on the identical layer accessed via PostGIS.

Also: LIKE works good with local shapefiles in Open Attribute Table > Advanced Filter (but this is not the solution - one can only select rows using a query but not filter a layer).
As we can read here, different behaviour between Query Builder and Attr. Table Advandced Filter comes from the fact that "one is handled by OGR, the other one by QgsExpression". If this is the case, then why does Query Builder work correctly on postgis layers?

I know this bug was reported several times (and even marked as 'wontfix'). Last time it was perhaps 4 years ago but nothing changed since then.

IU there any other way to filter shapefiles using SQL query?


texts: How to set the Identify Feature tool in QGIS 2.4?

texts:
The Identify Results window does not appear when using the identify feature tool  in QGIS 2.4.
I guess there's a check box somewhere I need to check/uncheck, but the settings of the tool changed compared to the previous version (see attached printscreen).


texts: What Datawarehouse can process geospatial informations?

texts:
I'm currently using a PostgreSQL database with the PostGIS extension. My querys are growing to big for my inhouse server, so I would like to port this into the cloud.

The querys mainly compute the distance between a table of geometrys and a table of points, returning all points within a certain distance to the geometry.

Does any cloud DWH (BigQuery, Redshift, ...) support processing geospatial information? If not, do you guys have an idea on how to go on with my problem?

edit: As requested by @JohnPowellakaBarça I'm posting the SQL query + explain that leads to my machine crashing.

explain
select t1.osm_id as building_osm_id, t2.osm_id as poi_osm_id,
t1.centroid::geography <-> t2.geom::geography as distance
from germany_buildings_centroid t1, germany_poi_classes_transformed t2
where t1.st_x <= (t2.st_x+0.01) AND t1.st_x >= (t2.st_x-0.01)
AND t1.st_y <= (t2.st_y+0.01) AND t1.st_y >= (t2.st_y-0.01)
AND ST_DWithin(t1.centroid::geography, t2.geom::geography,1000);

Gather  (cost=1000.00..48193265925353.62 rows=809798 width=24)
   Workers Planned: 2
   ->  Nested Loop  (cost=0.00..48193265722916.41 rows=337416 width=24)
         Join Filter: ((t1.st_x <= (t2.st_x + '0.01'::double precision)) AND (t1.st_x >= (t2.st_x - '0.01'::double precision)) AND (t1.st_y <= (t2.st_y + '0.01'::double precision)) AND (t1.st_y >= (t2.st_y - '0.01'::double precision)) AND ((t1.centroid)::geography && _st_expand((t2.geom)::geography, '1000'::double precision)) AND ((t2.geom)::geography && _st_expand((t1.centroid)::geography, '1000'::double precision)) AND _st_dwithin((t1.centroid)::geography, (t2.geom)::geography, '1000'::double precision, true))
         ->  Parallel Seq Scan on germany_buildings_centroid t1  (cost=0.00..1304445.57 rows=11660458 width=56)
         ->  Seq Scan on germany_poi_classes_transformed t2  (cost=0.00..265650.37 rows=7031637 width=56)

Indexes on my tables:

           table_name            |               index_name               | column_name
---------------------------------+----------------------------------------+-------------
 germany_buildings_centroid      | idx_german_building_centroid_centroid | centroid
 germany_buildings_centroid      | idx_german_building_centroid_geom     | geom
 germany_poi_classes_transformed | idx_germany_poi_classes_transformed    | geom

After changing the query:

  select t1.osm_id as building_osm_id, t2.osm_id as poi_osm_id
from germany_buildings_centroid t1, germany_poi_classes_transformed t2
where ST_DWithin(t1.centroid::geography, t2.geom::geography,1000);
                                                                                QUERY PLAN
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 Gather  (cost=1000.00..46609808612882.03 rows=65674991 width=16)
   Workers Planned: 2
   ->  Nested Loop  (cost=0.00..46609802044382.93 rows=27364580 width=16)
         Join Filter: (((t1.centroid)::geography && _st_expand((t2.geom)::geography, '1000'::double precision)) AND ((t2.geom)::geography && _st_expand((t1.centroid)::geography, '1000'::double precision)) AND _st_dwithin((t1.centroid)::geography, (t2.geom)::geography, '1000'::double precision, true))
         ->  Parallel Seq Scan on germany_buildings_centroid t1  (cost=0.00..1304522.12 rows=11668112 width=40)
         ->  Seq Scan on germany_poi_classes_transformed t2  (cost=0.00..265691.35 rows=7035735 width=40)


texts: Reloading aprx that has been modified using ArcPy with ArcGIS Pro

texts:
I'm writing a python script tool for ArcGIS Pro that uses a custom GUI created with pyqt. I would like the tool to be able to add and remove layers from the main map in the ArcGIS Pro project it is run from. This works fine normally, but I'm also trying to run the tool in a separate thread so that I can continue using other geoprocessing tools (see here). I can still add layers to the project and save the changes using my GUI, but since I am referencing the project by its full path (rather than arcpy.mp.ArcGISProject("CURRENT")), the new layers won't show up unless I close the project and open it again.
Is there anyway to 'reload' an aprx project that has been modified by another process without closing it first?


texts: Disabling/overriding internal arcmap commands with vba

texts:
I am creating a template for ArcGis JTX. This is a mxd file, with a vba onload event handler which removes some functionality from vanilla ArcMap.

I'm stumped on one problem: double clicking on a featurelayer opens a layer properties screen. I'd like to disable that functionality.

An ideal solution would be to replace the command itself with an ICommand that shows a messagebox with the text "Disabled by administrator"

Another (less than ideal) solution would be to remove the ICommand from the ICommandsCollection, but I haven't found such an object in the documentation.

As an alternative, I'd replace the doubleclick eventhandler from the featurelayer, disable shortcuts and remove buttons from toolbars and menus. It's a lot of programming, and I have yet to learn how to replace the eventhandler.

Do you have any idea if one of these solutions can be done?

tl;dr: I need a function DisableCommand(UID) in an mxd with vba.


texts: Time series chart of EVI (Landsat 8) for a single pixel in GEE

texts:
I'm trying to chart out EVI for a single pixel(roi) over time.
Basically trying to do what this link does for NDVI and apply it to the EVI available from landsat 8 imagery, but have no clue how to do it.
Heres my code for NDVI.
  var cliptocol = function(image){
  return image.clip(table);
};
var table = ee.FeatureCollection("users/napopovicc/nikanotee")

// Import the Landsat 8 TOA image collection.
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');

// Map a function over the Landsat 8 TOA collection to add an NDVI band.
var withNDVI = l8.map(function(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
});

// Create a chart.
var chart = ui.Chart.image.series({
  imageCollection: withNDVI.select('NDVI'),
  region: roi,
  reducer: ee.Reducer.first(),
  scale: 30
}).setOptions({title: 'NDVI over time'});

var cloudlessNDVI = l8.map(function(image) {
  // Get a cloud score in [0, 100].
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');

  // Create a mask of cloudy pixels from an arbitrary threshold.
  var mask = cloud.lte(10);

  // Compute NDVI.
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');

  // Return the masked image with an NDVI band.
  return image.addBands(ndvi).updateMask(mask);
});

print(ui.Chart.image.series({
  imageCollection: cloudlessNDVI.select('NDVI'),
  region: roi,
  reducer: ee.Reducer.first(),
  scale: 30
}).setOptions({title: 'Cloud-masked NDVI over time'}));


texts: How to get the geometry of all tables in a postgresql database?


texts: Molusce Value error: operands could not be broadcast together with shapes (13,)(14,)

texts: I am a PhD student and working on LULC changes. I am trying to use MOLUSCE. When I go to "Area Changes" tab, and I press the "Update tables" or "Create changes map", I always get the same python error message: -Python code: Value error: operands could not be broadcast together with shapes (13,)(14,). What can I do to solve this?


texts: SQL Server data types supported in ArcGIS

texts:
I am wondering if there is a known difference in processing speed in ArcGIS between the three geometry data types in SQL Server: compressed binary, SQL Server geometry, and SQL Server geography.

We have two SDEs in questioning on this issue

The first, and older, SDE has feature classes using Compressed Binary and showing the Geometry Property - "Storage: High Precision using Esri Binary spatial Type."
The second SDE, which was recently created, is almost an exact copy of the first SDE expect that the feature classes use SQL Server Geography and show the Geometry Property - "Storage: High Precision using Geometry spatial type (SRID 2278).

The only other difference's between the two SDE's is that the second SDE is in a new schema and domain tables were added.

Our users are reporting that the new SDE is running much slower than our original SDE we used to use while the feature classes are being edited.

Seeing as there could be many reasons as to why this is I figured I would start with the Geometry Properties because I did notice they were different between the two SDE's and as far as I know these two different storage types utilize spatial indexing in a different way which could impact performance.

Does anyone have any advice before we continue troubleshooting?


texts: How to connect to MySQL from QGIS in Ubuntu10?

texts:
I have just installed Quantum GIS 1.8.0 on my Ubuntu 10. I want to add a MySql table as a data source but I can't see any options for connecting to mysql.

Can any one guide me how to connect to mysql?

N.B : Options for Available Database Types in Qgis are MSSQL,ODBC,PostgreSQL and ESRI Personal GEO Database.


texts: GeoMedia + WFS 1.1.0 -- how to specify a non-default SRS in a GetFeature request?

texts:
I'm forming a WFS post request like the one below.  Notice I have SRSNAME="epsg:4326":

<GetFeature xmlns:wfs="http://www.opengis.net/wfs" service="WFS"
version="1.1.0" maxFeatures="500">
  <wfs:Query typeName="SewerMainline" SRSNAME="EPSG:3857">
    <Filter xmlns="http://www.opengis.net/ogc">
      <BBOX>
        <PropertyName>LinearGeometry</PropertyName>
        <Envelope srsName="EPSG:4326" xmlns="http://www.opengis.net/gml">
          <lowerCorner>46.226130203549 -119.23339223652</lowerCorner>
          <upperCorner>46.226323539414 -119.23301630567</upperCorner>
        </Envelope>
      </BBOX>
    </Filter>
  </wfs:Query>
</GetFeature>

But the result looks like this.  Notice the geometry is still in EPSG:3857 (the DefaultSRS):

<gmgml:FeatureCollection xmlns="http://www.intergraph.com/geomedia/gml"
xmlns:gml="http://www.opengis.net/gml"
xmlns:gmgml="http://www.intergraph.com/geomedia/gml"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation='http://www.intergraph.com/geomedia/gml http://usalvwdgis1:801/WFS1/request.aspx?request=DescribeFeatureType&amp;VERSION=1.1.0&amp;SERVICE=WFS&amp;OUTPUTFORMAT=text/xml;subtype=gml/3.1.1'
numberOfFeatures='4'>
  <gml:boundedBy>
    <gml:Envelope srsName="EPSG:4326">
      <gml:lowerCorner>46.2260082243399 -119.233511477885</gml:lowerCorner>
      <gml:upperCorner>46.2262617590092 -119.232662186944</gml:upperCorner>
    </gml:Envelope>
  </gml:boundedBy>
  <gml:featureMember>
    <gmgml:SewerMainline gml:id="SewerMainline.5635">
      <gmgml:ID1>5635</gmgml:ID1>
      <gmgml:LinearGeometry>
        <gml:LineString srsName="EPSG:4326">
          <gml:posList srsDimension="2">46.2262145488833 -119.233114496317
          46.2262139292097 -119.233327445878</gml:posList>
        </gml:LineString>
      </gmgml:LinearGeometry>
    </gmgml:SewerMainline>
  </gml:featureMember>
</gmgml:FeatureCollection>

I've confirmed the service supports 4326..here is a snippet from GetCapabilities:

<ows:Parameter name="srsName">
  <ows:Value>EPSG:3857</ows:Value>
  <ows:Value>EPSG:4326</ows:Value>
</ows:Parameter>

I also looked at the WFS 1.1.0 spec and the doc on the SRSNAME state the following:

  This parameter is used to specify a WFS-supported SRS that should be
  used for returned feature geometries. The value may be the DefaultSRS
  or any of the OtherSRS values that a WFS declares it supports in the
  capabilities document. The SRS may be indicated using EPSG codes or
  the URL form defined in [2]. If the parameter is not specified then
  the value of the DefaultSRS for the feature type being queried shall
  be used.

Where am I going wrong?


texts:
I have a need for a model that all the input may be optional and allow the "Selected features only" option.
So, I developed a model that fakes it all as seen in the example below

I realize that some SAGA algoritms acepts none input, so it can start de process allowing the selected features either.
I create some regular points to use in the merge, so the next algoritm will always receive something and don´t cause an error.

The question is.
Is this the only and viable solution, or the long way one?
because regardless the input, it takes at least 6 seconds to finish, even empty.


texts: Hive GIS ST_GeodesicLengthWGS84 not returning expected distance

texts:
I am trying to find the distance between two points (in meters and not degrees) using GIS Tools for Hadoop by ESRI. The ST_Length provides distance in degrees and hence I used the ST_GeodesicLengthWGS84 function. However, I noticed that the distances are not what I would expect to see. Below are the queries I tried:

SELECT  ST_GeodesicLengthWGS84(ST_SetSRID(ST_Linestring(0.0,0.0, 0.3,0.4), 4326))/1000;
SELECT ST_GeodesicLengthWGS84(ST_SetSRID(ST_Linestring(40.773773, -73.872307, 40.973773, -73.772307), 4326))/1000;
SELECT ST_GeodesicLengthWGS84(ST_SetSRID(ST_Linestring(40.773773, -73.872307, 33.944115, -118.40377), 4326))/1000;

I also used the distance calculator to check the distance (http://andrew.hedges.name/experiments/haversine/) and also google maps for directional reference (Google is a road distance, so it cannot be less than the distance from the above queries)

Here are the results I see for the three queries:
1) 55.42 (55.615 KM from distance calculator)
2) 12.77 (23.73 KM from distance calculator and 33.2 KM from google maps)
3) 4962.36 (3966.188 from distance calculator and 4536 KM from google maps)

Can someone tell me in case I misunderstood the use of this formula? And what would be the right way to get the distance in meters/KM?


texts: Using ArcGIS Field Calculator to change text to digit?

texts:
I am using ArcGIS Desktop 10.6.

I have one field that is titled SEVERITY, and is string type, and has words like "Medium" "Low" "High," etc.  I would like to assign a number to correspond with the word ("1" for "Low", 2 for "Medium," etc.), and to have that populate in a new field.

I've created the new field, called it "SEVERITY_CODE" and made it an integer.  All I need to do is use the field calculator to say "look at the SEVERITY column, and translate that into numbers for the "SEVERITY_CODE" column but I cannot seem to find it.  The solutions I've found on the web are not helping.

Can anyone help with this?


texts: How tow use an ArcGISCache layer along with an ArcGIS93Rest layer?

texts:
I want to have two layers in my OpenLayers map. The BaseLayer, as well as the overlay layer are being served by ArcGIS Server.

I have a problem with ArcGIS93Rest when using ArcGISCache as basemap.

The problem is ArcGIS93Rest doesn't show up and there is no HTTP request to the service.

Is it possible to use ArcGISCache and ArcGIS93Rest together?

Here is the code :

var baseLayer = new OpenLayers.Layer.ArcGISCache("AGSCache", "http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer", {
    isBaseLayer: true,
    resolutions: resolutions,
    tileSize: new OpenLayers.Size(layerInfo.tileInfo.cols, layerInfo.tileInfo.rows),
    tileOrigin: new OpenLayers.LonLat(layerInfo.tileInfo.origin.x, layerInfo.tileInfo.origin.y),
    maxExtent: layerMaxExtent,
    projection: 'EPSG:' + layerInfo.spatialReference.wkid
});
var forestLayer = new OpenLayers.Layer.ArcGIS93Rest("ForestLayer", "http://192.168.1.36:6080/arcgis/rest/services/forest_editing/MapServer/export", {
    layers: "0,1,2",
    TRANSPARENT: true
});


texts:
I am working with a number of detailed ArcGIS Geodatabases with domains used to enforce specific workflows, etc.  This is all good and great, but I would like to move my data into either a PostgreSQL/PostGIS or SQL Azure datastore and would like the ESRI client APIs (Mobile SDK, Web APIs, and Desktop templates) to enforce workflows in a similar manner to Geodatabase domains.  Obviously a domain is nothing more than a table constraint, but there is lots of built in ArcGIS functionality which makes life easier for the end user (mainly desktop users).  This seems very timely since it appears that at 10.1 ESRI is expanding support for native spatial support (i.e. Simple Features).

Will ESRI provide domain support within a native spatial DB from ArcGIS Desktop?  Thoughts, suggestions, etc.


texts: Is it possible to import a Shapefile from Manifold into QGIS with the same projection?

texts: is it possible to import a shapefiles from manifold to Qgis with the same projection?


texts: How to check if a point is snapped to a line in QGIS

texts:
I have two vector classes, sg_point and sg_line. Any sg_point HAS to be snapped to a sg_line. I should not find errors if I use snapping correctly, but in case there are errors how do I check in QGIS?

I have tried select by locaction, but no matter what option I use, either I get all points (useless) or no points (also useless). I cannot get either only the points that are snapped nor only the points that are not snapped to a sg_line instance.

What am I doing wrong?

I am using Pisa 2.10.1

UPDATE -
I have tested many configurations. It seems to me that select by location only works when the snapped point is not a vertice. When it is snapped to a segment, it is not selected. Is this working as expected, or is this a bug?


texts:
I have 3 SpatiaLite layers (one without geometry but with fields that need to be filtered, one with points and one with polygons, both without any non-geometry fields) that need to be combined to a single one, i.e. a left outer join is needed. Also, some filtering on the attributes needs to be done (in particular, a certain time range needs to be filtered). Apart from that, I also need to spatially filter the layers with points and polygons:

All the points inside a certain boundary (a single polygon only available as a SHP) need to be retained
All the polygons that overlap with this certain boundary need to be retained, i.e. both the polygons fully inside the boundary and also those which are only partially inside.

The main problem is that all of these layers contain quite a lot of data (i.e. >800mb) and it is quite tedious to trial-and-error the different possibilities. That is why I am asking what work process I should follow:

Filtering and joining the SpatiaLite data with the Database manager in QGIS and then doing the clipping or whatever geometry function is needed for the above process? Could I actually do the clipping or whatever geometry function also directly in SpatiaLite, given that the boundary is only available as a SHP?
Doing the clipping or whatever geometry function is needed first, which yields a shapefile, but how do I join them afterwards?

Actually, I am pretty lost / don't know the QGIS functions well enough.


texts: Is there a way to select all items within a circular buffer?

texts:
I want to select all the items within a circular buffer. Inside the buffer are
multipolygons and multipoints. See the screenshot i have posted.


texts: Recommended way for changing value in PostGIS using GeoServer

texts:
I have a map being rendered from PostGIS using GeoServer. I frequently make some edits to the map and update the new value by directly connecting to the PostGIS database and then making the changes using SQL query.

Is this the recommended way to do this or does GeoServer offer any other proper way of achieving the same ?

I read slightly about WFS, WCS and WMS and doubt whether they are designed for this task, as I always get lost in their documentation.


texts: Returning a multipart polygon back to singlepart after editing

texts: Feature "A" had two polygons that would show up together as one record but after editing it turned into two separate records. I would like to get them back together so they are identified together and have one single record. What can I do to make that happen?


texts: How to add color interpretation for raster bands using GDAL?

texts:
I have 3 separate 1-band GeoTIFF files. The color interpretation for band is gray. I want a 3-band RGB file. I've used gdal_merge.py

gdal_merge.py -separate file1.tif file2.tif file3.tif -o output_file.tif

But the color interpretation for the 3 bands of the resulting output_file.tif is gray, undefined, undefined. Apart from that, all pixels are 0.

If I execute pct2rgb.py, I get a RGB file:

pct2rgb.py output_file.tif output_file_rgb.tif

But of course, the pixels are 0 too. So, I have 3 questions:

Is gdal_merge.py the right tool to combine 3 1-band files in one 3-band RGB file?
Why am I getting undefined color interpretation for bands?
Is pct2rgb the right tool to transform 3-band files with this color interpretation to RGB files?

UPDATE: The rasters does not have a color table. Just color interpretation: Gray.

On other hand, pixel values goes from 0 to 1023 (this is deliberate)

More data: they're rotated rasters (no north up), but all of them have the same geotransform.

UPDATE 2: I can warp the images to make them north up, construct a VRT and add ColorInterp for each band, but I still get color interpretation as gray, undefined, undefined in the output result.

The problem is I need to create a color table in, at least, the first band. I know the way to create them, but I don't know how many entries should my table have. Why are there 13 entries in the example of GDAL Raster FAQ? All the pixels have values between 0 and 1023, if helps.

UPDATE 3: Apparently, there's no way in the TIFF format to really specify the color interpretation of each band. The way GDAL builds the color interpretation when reading a TIFF file is a combination of the value of the PHOTOMETRIC and EXTRASAMPLES tag.

Reading about these tags:

PHOTOMETRIC represents the color space of the image data. A value of 2 means that the components of a pixel value are RGB, but it assumes Byte pixels, and I have UInt16 pixels (I tried -co "PHOTOMETRIC=rgb", and got an error). So, I can't specify PHOTOMETRIC tag for the output file.
EXTRASAMPLES specifies that each pixel has N extra components. I'm not sure about how to use this tag to create my merged file. Or if a I need it.

So, in update 2 I suggest the creation of a ColorTable, but how? In my 3 input files, pixel values go from 0 to 1023. Do I have to match them with colors? Do I have to create a ColorTable with 1024 inputs? How?

In update 3, seems that I could use some GeoTIFF tags when creating the merged file, but I'm not sure if I really can use them, or how.


texts: Interpretation of historical aerial imagery

texts:
I've been having to look at historical aerial imagery with my job as of late. The goal is to identify and describe historic land use of the sites. Specifically, I've come across sites that were previous ag fields that have a common element that I've been unable to identify.

It looks like little dots in concentrated areas of the landscape, possibly trees, maybe stumps. In future photos, they are often not present suggesting that they are not living trees at least.

I've attached a photo with a couple example sites outlines by a green line.

Any thoughts?


texts: Plotting classification map using shapefile in matplotlib?

texts:
"I have been struggling with python map plotting since past few days. I have realized that one can use python for plotting comparatively easy maps, but when one has to display detailed information, it doesn't work."

I am plotting land use land cover map in python. I have 8 different classes represented using 1,2,3,4,5,6,7,8. Each class need to be shown with unique color.
I am not getting the 8 different colors for 8 different classes, as it takes the number while creating colormap. I have attached my map below

Here the classes are 8, but the color codes are 7.

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
from mpl_toolkits.basemap import Basemap

num_colors = 8
cm = plt.get_cmap('jet')
scheme = [cm(i*1./num_colors) for i in range(num_colors)]
bins = np.linspace(1, 8, num_colors)

fig = plt.figure(figsize=(18, 10))
ax = fig.add_subplot(111, axisbg='w', frame_on=True)

m = Basemap(llcrnrlon=85.9,llcrnrlat=22.,urcrnrlon=87.9,urcrnrlat=23.6,
             resolution='i', projection='tmerc', lat_0 = 22., lon_0 = 87.)

parallels = np.arange(15.,30.,0.25)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=12, linewidth=0.5)

meridians = np.arange(70.,90.,0.25)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=12, linewidth=0.5)

m.readshapefile('Path\To\Shapefile', 'units', color='#444444', linewidth=0.5)

for info, shape in zip(m.units_info, m.units):
    LULC = int(info['LULC_Recod'])
    which_bin = np.digitize([LULC], bins)
    color = scheme[which_bin-1]
    patches = [Polygon(np.array(shape), True)]
    pc = PatchCollection(patches)
    pc.set_facecolor(color)
    pc.set_edgecolor('none')
    ax.add_collection(pc)

ax_legend = fig.add_axes([0.8, 0.1, 0.01, 0.8])
cmap = mpl.colors.ListedColormap(scheme)
cb = mpl.colorbar.ColorbarBase(ax_legend, cmap=cmap, ticks=bins, boundaries=bins, orientation='vertical')
cb.set_ticks([i for i in range(1, 9, 1)])

cb.set_ticklabels(['Water', 'Deciduous Broadleaf Forest', 'Mixed Forest', 'Woody Savannas',
    'Grasslands', 'Croplands', 'Urban and Built-Up', 'Cropland/natural vegetation mosaic'])

plt.show()

How to plot this classification map correctly?


texts: Apply a symbology from a style to a layer in a map document using arcpy?

texts:
I've looked at a number of sources and from what I gather, it seems like the only way I can apply a symbology is to do the following:

Manually create a layer file that has the symbology I want already applied.
Apply the symbology from the created feature layer using the arcpy.ApplySymbologyFromLayer_Management function.

This is what I have in code so far to do the above:

mxd = arcpy.mapping.MapDocument(MAP_PATH)

data_frame = arcpy.mapping.ListDataFrames(mxd)[0]

for layer in arcpy.mapping.ListLayers(mxd):
    symbol_name = path_to_symbol_layer_file
    temp_layer_name = layer.name

    arcpy.MakeFeatureLayer_management(
        in_features=layer.dataSource,
        out_layer=temp_layer_name,
    )

    arcpy.ApplySymbologyFromLayer_management(
        in_layer=temp_layer_name,
        in_symbology_layer=symbol_name,
    )

    arcpy.mapping.UpdateLayer(
        data_frame=data_frame,
        update_layer=layer,
        source_layer=arcpy.mapping.Layer(temp_layer_name),
        symbology_only=True,
    )

mxd.save()

del mxd

I'd prefer to keep my symbologies organized in a single .style file and apply each one directly, without having to create a layer file for it. Is this not possible with arcpy?

In other words, in pseudo-code, I want something like:

style = arcpy.Style('/path/to/foo.style')

mxd = arcpy.mapping.MapDocument('/path/to/bar.mxd')

for layer in arcpy.ListLayers(mxd):
  symbology = arcpy.GetSymbology(style, "custom_symbology_%s" % layer.name)

  arcpy.ApplySymbology(layer, symbology)

mxd.Save()

del mxd


texts: Modify interaction using Icon image in OpenLayers3

texts:
I have vector features represented by icon images in a map using OpenLayers 3 (3.0.0). These features can be modified by the Modify interaction only if the user place the mouse cursor close their position. If the user tries to drag the icon, outside the position range, the map is panned (instead of dragging that feature).

Is there any way to expand the drag area, from the Modify interaction, to the feature icon?

Example: http://jsfiddle.net/suguiura/xuzsj99b/


texts: DTM to TIF iteration asking for projection

texts:
I have generated a script where all the .dtm files in a directory will be translated into .tif in another directory.

The issue here is that in some computers, a window appears after each iteration asking for projection (even if the tif file is not added to the canvas) and the process stops until you close that window.

I guess it must be something to do in order to avoid this window appearance, almost during an script execution.

I use to work with qgis 2.18.12


texts: Can I import Network Dataset in PostgreSQL through pgRouting or Polyline Shapefile?

texts: I need to make an application for the Quickest Route to the nearby hospitals from the incident place. I have read many articles and papers but I am confused whether I need polyline shapefile or network dataset or any other format for storing data in PostgreSQL


texts: Find intersecting polygon within a layer

texts:
I would like to find all polygons, which intersect another polygon within the same layer. I thought the geometry checker could do the job, but it doesn't.

All polygons are intersecting except one, but I get no geometry errors. What I am doing wrong? (QGIS 3.8)


texts: JTS GeometricShapeFactory generate an ellipse properly

texts:
This is what use to get a circle:

GeometricShapeFactory shape = new GeometricShapeFactory(gf);
        shape.setCentre(new Coordinate(40.748754,-73.985460));//Empire Sate Building
        shape.setSize(2* (10*0.009));//expect 10km radius
        shape.setNumPoints(64);
        final Polygon circle = shape.createEllipse();//createCircle deprecated

What I get is a weird ellipse. I know that closer to the equator we get 'more perfect' circle but still, is it valid circle for radius of 10km with ESB as a center point ?

Reading Dr. Ian Turton posts, github repo, and other great materials he shares, I know that I can use geotools and do that differently but I would like to stick with that approach.

EDITED:

I am writting an application. In MySQL I store malls with their location (spatial type Point). What I want to do is to get malls from db that are in a certain circle (via ORM). In order to do that I have to create a Polygon(circle, linearring) in which orm will look. And Instead of using complicated calculations (and additional library). I thought I can use GeometricShapeFactory from jts lib I am already using. But because I am less than beginner what it comes to these things, all I wanted to ask, is whether or not the vertical ellipse generated by jts is valid for my purposes or not ?

I am using google maps (as basemap).

My question:
Based on the fact that my ellipse is just wrong (considering given radius it should be a perfect circle - @whuber comment), how can I generate a proper Polygon(circle, linearring) based on radius and center point ? (I am going to use it in hibernate-spatial: (...).add(SpatialRestrictions.within(Mall.LOCATION, getCircle(centerPoint, 10)));)


texts: Zero-size array to reduction operation maximum which has no identity

texts:
Below is my identify background code, which identifies a background from an image and removes it.;

Previously it's was working fine but after I chunked my images and try to find out background I got an error.

  Traceback (most recent call last):

File "create_vips_test.py", line 336, in

print test_h5_file()

File "create_vips_test.py", line 243, in test_h5_file

res =  identify_background(tile).astype(np.int16)

File "create_vips_test.py", line 111, in identify_background

max_area=np.max([cv2.contourArea(cnt) for cnt in contour])

File "/usr/local/lib/python2.7/site-packages/numpy/core/fromnumeric.py",

line 2293, in amax
out=out, **kwargs)

File "/usr/local/lib/python2.7/site-packages/numpy/core/_methods.py",

line 26, in _amax
    return umr_maximum(a, axis, None, out, keepdims)

ValueError: zero-size array to reduction operation maximum which has no identity

def identify_background(img):
im  = cv2.cvtColor(np.array(img), cv2.COLOR_RGBA2GRAY)
# print 'identiyf background ', im
blur = cv2.GaussianBlur(im,(5,5),0)
# print 'blur', blur
# find normalized_histogram, and its cumulative distribution function
hist = cv2.calcHist([blur],[0],None,[256],[0,256])
hist_norm = hist.ravel()/hist.max()
Q = hist_norm.cumsum()

bins = np.arange(256)

# print 'bins', bins
fn_min = np.inf
thresh = -1
for i in xrange(1,256):
    p1,p2 = np.hsplit(hist_norm,[i]) # probabilities
    q1,q2 = Q[i],Q[255]-Q[i] # cum sum of classes
    b1,b2 = np.hsplit(bins,[i]) # weights

    # finding means and variances
    m1,m2 = np.sum(p1*b1)/q1, np.sum(p2*b2)/q2
    v1,v2 = np.sum(((b1-m1)**2)*p1)/q1,np.sum(((b2-m2)**2)*p2)/q2

    # calculates the minimization function
    fn = v1*q1 + v2*q2
    if fn < fn_min:
        fn_min = fn
        thresh = i

# find otsu's threshold value with OpenCV function
ret, otsu = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
# print 'ret', ret
# print 'otsu',otsu
orig=otsu
final=cv2.bitwise_and(orig, cv2.bitwise_not(orig))
# print 'final', final
im_bw_inv = cv2.bitwise_not(orig)
# print 'im_bw_inv', im_bw_inv

contour, _ = cv2.findContours(im_bw_inv, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)
# print 'contour, _', len(contour), _
max_area=np.max([cv2.contourArea(cnt) for cnt in contour])
# print 'contour', contour
# print 'max area', max_area

for cnt in contour:
    cv2.drawContours(final, [cnt], 0, 255, -1)
# print 'final', final
return final


texts:
Anyone know where or how to obtain a database of worldwide beaches?

Required fields are:

name
latitude
longitude
coastline aspect (if possible)

Alternatively, is there a way to query Google Maps to find the closest beaches to a given geolocation?


texts: Geoserver does not wrap GRIB raster points to coordinates where x < 0

texts:
I am new to the GeoServer and want to use it for plotting wind barbs based on a GRIB file over land data.

The wind barbs are properly created but nothing is shown on coordinates x < 0.
The map data has bounds (for x) of -180 - 180 whereas the grib data contains values between 0 - 360. Shouldn't the Geoserver be able to handle the wrapping? At least for the land data it works...

Here is an image of the results I get:

The strange thing is, that when I zoom out completely, the wrapping seems to work fine:

Can you hint me in any direction whether my data is wrong, the layer is not defined correctly or if it is an issue with the style declaration?

EDIT:

In my style definition for the wind barbs, I specify the spatial resolution for different zoom levels for the raster  as follows:

<ogc:Function name="parameter">
    <ogc:Literal>scale</ogc:Literal>
    <ogc:Function name="Categorize">
        <ogc:Function name="env">
            <ogc:Literal>wms_scale_denominator</ogc:Literal>
        </ogc:Function>
        <ogc:Literal>64</ogc:Literal>
        <ogc:Literal>100000</ogc:Literal>

        <ogc:Literal>32</ogc:Literal>
        <ogc:Literal>500000</ogc:Literal>

        <ogc:Literal>8</ogc:Literal>
        <ogc:Literal>1000000</ogc:Literal>

        <ogc:Literal>2.0</ogc:Literal>
        <ogc:Literal>5000000</ogc:Literal>

        <ogc:Literal>0.8</ogc:Literal>
        <ogc:Literal>10000000</ogc:Literal>

        <ogc:Literal>0.4</ogc:Literal>
        <ogc:Literal>20000000</ogc:Literal>

        <ogc:Literal>0.3</ogc:Literal>
        <ogc:Literal>60000000</ogc:Literal>

        <ogc:Literal>0.1</ogc:Literal>
        <ogc:Literal>70000000</ogc:Literal>

        <ogc:Literal>0.07</ogc:Literal>
        <ogc:Literal>140000000</ogc:Literal>

        <ogc:Literal>0.05</ogc:Literal>
    </ogc:Function>
</ogc:Function>

Also I found out, that when I convert the grib file into a geoTiff with

gdalwarp -t_srs EPSG:4326  inputGrib.grb output.tif --config CENTER_LONG 0 -wo SOURCE_EXTRA=100

everything works as expected.

EDIT 2020-02-26
This behavior is caused by a bug as described on the geotools issue tracker.
It will be fixed in GeoTools v22.3 which is used in Geoserver v1.17 will be released in March 2020.
Unfortunately, wrapping still fails for wms requests with output format application/pdf and image/svg+xml.


texts: Should I rasterio.open the destination when reprojecting raster?

texts:
I use the code:

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

mask=rasterio.open(path1)
dst_crs='EPSG:4326'
with rasterio.open(path2) as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, mask.width, mask.height, *mask.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
        })

    with rasterio.open(pathdst, 'w', **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

that gives now an error RasterioIOError: Blacklisted: file cannot be opened by driver 'netCDF' in 'w' mode.

My src file is multibanded, so I wat to reproject all bands at a time. Mask file should give all information for transformations. What is the shortest and the most pythonic way to do that?
Note: my files have no SRS (it is None).

After comment:

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from netCDF4 import Dataset

mask=rasterio.open(path1)
dst_crs='EPSG:4326'
with rasterio.open(path2) as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, mask.width, mask.height, *mask.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
        })

    dst=Dataset(pathdst, "w", format="NETCDF4")
    for i in range(1, src.count + 1):
        reproject(
            source=rasterio.band(src, i),
            destination=rasterio.band(dst, i),
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=transform,
            dst_crs=dst_crs,
            resampling=Resampling.nearest)

so what is the way to adress to bands with netCDF4? There is nothing about it in the docs.
Of course it gives an error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-102-dbd35e15c61a> in <module>
     21         reproject(
     22             source=rasterio.band(src, i),
---> 23             destination=rasterio.band(dst, i),
     24             src_transform=src.transform,
     25             src_crs=src.crs,

/opt/conda/lib/python3.7/site-packages/rasterio/__init__.py in band(ds, bidx)
    251     rasterio.Band
    252     """
--> 253     return Band(ds, bidx, set(ds.dtypes).pop(), ds.shape)
    254
    255

netCDF4/_netCDF4.pyx in netCDF4._netCDF4.Dataset.__getattr__()

netCDF4/_netCDF4.pyx in netCDF4._netCDF4.Dataset.getncattr()

netCDF4/_netCDF4.pyx in netCDF4._netCDF4._get_att()

netCDF4/_netCDF4.pyx in netCDF4._netCDF4._ensure_nc_success()

AttributeError: NetCDF: Attribute not found


texts: How to compare two different polygon shapefiles to compare how shapes have changed between the two?

texts:
I am trying to use R to compare the changes in shape of two multipolygon shapefiles. Specifically, these data are congressional districts from 2010 and 2013 in Wisconsin, so they have the same number of rows, with different shapes.

I am curious how I would compare these to each other to see how they have changed, similar to this question previously asked on here: Identifying differences in polygons across shape files.

So far I have only been able to read both shapefiles into R using the sf and rgdal packages.


texts: Fitting to Margins of Layout in ArcObjects?

texts:
What ArcObjects are involved in making "Fit to Margins" work for a given Page Layout?

Easier to see from picture below:


texts: Creating line with three points with Python in QGIS?

texts: I'm new to Python and I'm having some dificulties. I want to create a simple layer with a line in the console of QGIS with Python. How can I do that?


texts: How can I merge two geometry columns of different tables to a single column in another table

texts:
I have many tables with geometry attribute which is saved in Postgis database. These tables are basically from more than one shapefiles which is stored in Postgis using geotools.

Now I want a common table to handle all the geom attributes in the table, as shown in the below example

For eg:

Is there a way to do this using PostGIS? Initially the Table geom_all with geom column contains NULL value. I tried updating the values using JOINS with the below query example.

 UPDATE Geom_all g SET geom = (SELECT geom FROM Table-1 t WHERE g.id = t.id);

But this overrides with every new value i update in the geom column and also the sub query should return only 1 row, which in my case is not. I am not sure whether I can achieve this task using JOINS, it would be helpful to know if there is any other way to solve this.

NOTE: All the geom column contains geometry_type as MULTIPOLYGON and ST_SRID as 4326.


texts:
I have a toolbox to share with my collegaues, it includes a model and imported python script. This script is also imported to the model in the end of session.
When I share this toolbox to another PC, there is a problem in it: the model seems to be broken like here:

The solution is to edit model and manually set the path to a script, which is in the same toolbox or drag and drop it to a model.
Is there any kind of tool or calculator that will run my script without any relationships to location of script?


texts: Where can I download a raster dataset with the latitude and longitude of each pixel in geographic projection?

texts:
I want to do raster calculations in QGIS based on the latitude and longitude of each pixel.  The calculations will compute the area of the pixel, but also do other calculations for various projections of land masses.

The QGIS raster calculation facility does not appear to support variables for the X and Y coordinates (if it did, I could simply calculate the lat and long from the raster file metadata).

My raster maps are at 0.5 minute resolution, so a grid dataset of this resolution or better would be desirable.

The existence of many related questions on StackExchange indicate to me that there is a general need for a file such as I describe.  On the other hand, maybe there is an easy way to create such a file within QGIS.  That would be fine too.


texts: Leaflet Custom TileLayer - how to convert real-world coords to my map where centre is 0,0

texts:
I have a working leaflet map which utilizes a single TileLayer to display my custom map tiles.

It's a cartoon style map of a real location and everything is to scale. In fact you could place a google satelite map over the top of it and everything would line up perfectly.

I want to be able to plot real life lat/lon coords on my cartoon map.

The problem is, I've discovered it adopts it's own scale and coordinate system where the centre is 0,0.

I'm wondering if there's the possibility to do one of the following:

a) Specify in leaflet what the real life centrepoint and bounds of my map are?

OR

b) From marrying up various coordinates e.g: centrepoint and bounds of my cartoon map to the real life equivalents, devise some kind of formula to convert between the two.

I've actually spent some time on approach (b) but not been successful.
I thought I could do something along the following lines:

    latlngToXY: function(lat, lng) {

        // 0,0 on the cartoon map, marries up to these coords
        var mapCentre = {
           lat: -23.9076758,
           lng: 153.3127296
        };

        // Subtract mapCentre to get offset from 0
        var x = lat-mapCentre.lat;
        var y = lng-mapCentre.lng;

        // Get the scale multiplier, calculated from (lat-mapCentre.lat) divided by an already known cartoon map lat
        var scaleRatio = -30652.911561212319766;

        x = x*scaleRatio;
        y = y*scaleRatio;

        return [x,y];
  }

It worked fine for one of my coords (unsurprisingly the one I calculated the scaleRatio from), but is wildly inaccurate for anything else.

Maybe someone could point out what I'm missing or perhaps indicate if it's even going to be possible.


texts: Which hardware configuration is necessary for QGIS?

texts:
What kind of hardware is necessary to run QGIS on Win XP?
Closely related:
Specifications for GIS desktop hardware


texts: How to compute Travelling Salesman Problem by ArcGIS?

texts: I'm newbie in geography. I have a number of x-y coordinates in a text file. I want to compute the shortest path going through these points in ArcGIS. I can load my data points in ArcGIS, but network analyst options are disabled. Would you please help me how to do the computation?


texts: Calculation of fractional cover for each vegetation class at 30-m resolution matching MODIS swath in R

texts:
Currently I have a Landsat land cover map with 30-m resolution and I want to calculate the fractional cover for each vegetation class (forests, grasslands and croplands etc.) matching MODIS swath and do a linear regression analysis linking land cover percentage with the corresponding MODIS LST for each 1-km resolution pixel. Therefore, the fractional changes will be calculated by computing at 30m the proportion of each vegetation class (e.g. forest) that occupy each MODIS 1km grid cell.

I hardly found any practical tutorial using R to achieve this.


texts: How to publish python toolbox with lots of tools in ArcGIS Server 10.3?

texts:
I have python toolbox with like 40 tools and I need a quick way to publish the whole toolbox. I found some instructions in the documentation, they state that every tool should be executed and their results should be published by using ArcMap/ArcCatalog, but this will take forever!!! Because some tools accept maps in JSON, have a lot of complex parameters, needs extremely accurate parameters to execute succesfully, etc!

Is there a new simple way to publish the whole toolbox?


texts:
I've compiled and installed DARP branch. Is there any data example to test it?.

I've used some data based on document instructions http://www.pgrouting.org/docs/1.x/darp.html but no results return.

What's the difference between depot_id and depot_point_is?
In distance query i don't understand fields from_order and to_order.
If i add a request for going from point A to point B, is necesary to add in distances table the cost of going from point B to point A?.
The field values in distances query, what time units should be?, minutes?, seconds?, miliseconds?.
Where do darp.c and darp_solver.c leave log messages?.

I know that maybe darp branch is a RC version, but please could you help me with this?.


texts: Extracting the contents of a Google-KML file

texts:
I finally found all the datasets I could possibly need from http://globe.six.nsw.gov.au/ ...but a few hours later I'm not sure if it's possible to access this data into QGIS.

I've tried MapWindow (w/KML2Shapefile), GPSBabel, GlobalMapper and some online tools (as mentioned in other posts); I realise the file is essentially empty and contains a link but I can't seem to 'explore' that link any further.


texts: I am wondering if there are open source applications that make extensive uses of the special functions (the ST_* stored procedures) provided by PostGIS? I found many apps that support storing map data in PostGIS, but they seem to just fetch and store data from PostGIS and then operate on them in application code. I would like to learn about how to use those special functions as part of the queries by looking at sample applications. Thanks!


texts: How map unit affects the caption size in the label in QGIS?

texts:
I tried to make a the label caption size adjusted relatively to the map scale.
It worked, more or less, when I set the size into map unit.
I set them to 0.1, and when I zoom out the caption size get smaller

And when I zoom in it get bigger

However I cannot understand what is the "map point" logic, or in other words how does it work.
And one more question: is there a better way to adjusted the caption size and force them to always be inside the polygon?


texts: GeoDjango -How to store measurement data

texts:
Assume I have recorded the following values during a bike trip. One such dataset has been recorded every 5 seconds:

timestamp, longitute, latitude, pulse rate, pedal frequency

What is the best way to store them in a GeoDjango application? I thought about using a model similar to the following.

class Measurement(models.Model):
    comment = models.CharField()

class MeasurementPoint(models.Model):
    measurement = models.ForeignKey(Measurement)
    timestamp = models.DateTimeField()
    point = models.PointField(srid=4326)
    pulse = models.IntegerField()
    pedal = models.IntegerField()

This works and I can query the measurement points. But now I want to display the bike trip on a map, so I have to convert the individual points to a LineString. How is this done? And how would it bit possible to use different colors to highlight for example the pulse rate with different colors?


texts: How to create a map-based survey with drag-drop capability?

texts:
I am an urban planning master student. I need to create a map-based survey for my research. The survey should allow users to drag some different points from the toolbar to the map showing next to it. When a user finishes the survey coordinate data and types of the dragged points should be store. To make it clear i would like to collect data of the parks which people think it is safe, fun, better etc. To allow people to rate these parks and collect geolocated data using point generation seems logical. Here is a good example for this kind of survey, http://mydowntown.sasaki.com/houston/survey.html?relation=live . Using map API's like MapBox, Google Maps, or libraries like leaflet etc. is the most suitable way i suppose. I also have access to ArcGIS products, which are quite efficient when it comes to surveying. However my coding knowledge is like zero and im having hard time while scripting this kind of interaction without any basic example tutorial.

So how do i create a survey like above? Is there any pre-made applications for it, or ant examples that i can tweak?


texts: ArcMap Exported (KML/KMZ) Meridians and Parallels do not match/ align with that of Google Earth

texts:
I was trying to draw kml that are (one degree in gap) meredians and parallels for Austrailia. e.g see i tried to draw 142 degree and 11 degree lines using fishnet in arcmap but when imported into Google Earth they shifted a bit. I used WGS84(EPSG:4326) when creating lines.I supposed that google earth uses EPSG:4326 (I also tried EPSG:3857) as described at here, here and here.

Now how can i achieve 100% matched kml in arcmap. My extent is 100E, 50S to 155E, 3S.

My KML

Below are the steps i did to generate KML as attached above-

Open Create Fishnet Tool in a dataframe set to WGS84(EPSG:4326) and assign the
parameters as -

My fishnet setting-

From the spatial extent i am generating fishnet on i.e. 100E, 3S to 155E, 50S i counted the number of parallels(50-3=47) and meridians(155-100=55) and used this number as row and column number in the fishnet tool as shown above.
Then i used Densify tool to generate sufficient points in the lines(0.001 degrees distance interval) whereas projection does not shifts the whole feature.
I defined the projection of the layer WGS84(EPSG:4326) using Define Projection tool(I also tried EPSG:3857)
Finally i reprojected the line layer into WGS84(EPSG:4326) using Project tool(I also tried EPSG:3857)
And exported the lines into kml using Layer to KML tool.

Now, Problem is when i zoom in the exported kml in google earth i see a shift as shown in the attached.So i need a kml(with generation method) in the above spatial extent(100E, 50S to 155E, 3S) that matches exactly with the grid of google earth irrespective of zoom level

System datails:

I am in Windows 8 and Arcmap Advanced 10.3.0.4284 and by WGS84 i mean as below

GCS_WGS_1984
WKID: 4326 Authority: EPSG

Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
  Spheroid: WGS_1984
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314245179
    Inverse Flattening: 298.257223563


texts: Overpass API query that returns amenities and filters those without a name tag

texts:
I'm trying to query all the amenities in my neighborhood but find that the result set has a bunch of amenity-nodes without a name tag.

I want to filter these nodes out.

After reading the Overpass Language reference wiki (here) - http://wiki.openstreetmap.org/wiki/Overpass_API/Language_Guide#Tag_request_clauses_.28or_.22tag_filters.22.29

where it says

["key"!~"value"] filter objects tagged with this key but a value not
                 matching a regular expression

The following returns nothing.

/*
This has been generated by the overpass-turbo wizard.
The original search was:
“amenity=*”
*/
[out:json][timeout:25];
// gather results
(
  // query part for: “amenity=*”
  node["amenity"];node["name"!~"^$"]({{bbox}});
  way["amenity"]({{bbox}});
  relation["amenity"]({{bbox}});

);
// print results
out body;
>;
out skel qt;


texts: UTM coordinates from fishnet label points not correct in Google Earth UTM

texts:
I am a wildlife researcher and I am trying to create a set of points 1 km apart on a map of a particular area of my study region. I created a polygon (in Google Earth) that is about the size of what I need my study site to be over the study area, imported this into ArcMap and created a fishnet to the extent of the polygon (with the grid cells 1x1km).

With the fishnet function, evenly spaced points already come up, so what I want to do with those is extract the xy coordinates so I can put them into my GPS device and find them on the ground at my study site!

So, I followed instructions and extract values to each point and used my satellite image data layer of my map with my study site as the input raster; the geographic coordinate systems for all of my files are GCS_WGS_1984 and the projected coordinate system is WGS_1984_Costa_Rica_TM_90 (as my study site is in Costa Rica). I then did the add XY coordinates tool and entered the points shapefile I just created with the extracted raster values in the attribute table. However, the UTM coordinates I am given are way off my study site when I put them into Google Earth to check. I'm not really sure what is going on. I'm wondering if it is off because it is derived from the satellite image in the projected coordinate system rather than the geographic? Maybe it's skewed? Does anyone know how I might be able to fix this mismatch?

I would just use Google Earth but the point function in ArcMap is critical to spacing out my sampling sites correctly. I decided to test some other GPS coordinates I had added to ArcMap as XY data and they too were very off in Google Earth when put into UTM in ArcMap.

For example, one point in ArcMap was: 412327.28387 M E 1136639.13284 M N
and the same point in Google Earth was: 741099.72 m E 1137197.38 m N

However, one other thing is I do not know what the UTM zone is; Google Earth said 16 P, but ArcMap does not specify.


texts: Geoserver WMS image points drift to inaccurate position as I zoom out

texts:
I have an application that uses Leaflet maps and I've overlayed some geoserver WMS
layers (both point and Polygon). Now when I'm in at zoom 8 or more, etc
the overlay looks accurate and exactly where it should be - but as I move
out to zooms 7 and less I find my WMS layer seems to shift up more and more with each zoom level I go out.

The projection of the data is 4326 (WGS84). No errors, but as I said as you zoom out you notice the
gradual shifting.

I am attaching an image for clarity. Can anyone help me solve this problem?

---------------------Update-------------------

Please find the url below. This url returns an image that gets overlayed on the map. It is at zoom level 4. Its a very elaborate code & I could not put it here. Hence, I am giving the link that gets generated through the code and the returned image gets overlayed on the map. Every time a pan or zoom happens a new image as per the new view port gets requested and the old one gets replaced by the newly requested one.

We have been using our own proprietary layer imagery generation code in the past. Our data used to be stored in SQL Server as geography data type. We never faced this problem. In case of boundless geoserver we can only work with geometry data type (geography is not supported). Now when we are substituting our proprietary layers with boundless layers we are seeing this problem. Not sure if there is a setting or something that I need to do while setting up the layers in boundless geoserver that would take care of this issue.

http://aaaa.bbbb.cccc.com:8080/geoserver/usa/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&STYLES&SRS=EPSG%3A4326&LAYERS=usa%3ALayer1&BBOX=-152.40234375%2C27.839076094777816%2C-72.158203125%2C61.73152565113397&WIDTH=913&HEIGHT=569


texts: Setting current group as active group in PyQGIS?

texts:
While iterating on groups, I would like to set the current group as “an active group” in order that generated layers will be children of this active group.

Is there any pyqgis expression to do this (like for layer: iface.setActiveLayer(vl))?


texts: I have downloaded the Time Manager plugin for QGIS.  I have created a simple test file with only 4 points, with a date attribute showing YYYY-MM-DD, each point moving forward one day.  Yet it tells me that the attribute selected for use as start time contains invalid data. What am I doing wrong?


texts: Arcpy connecting points with lines using coordinates

texts:
Based on another post, I wrote a Python script that reads coordinates, creates points, and then connects them. Here is the code:

import arcpy

in_rows = arcpy.SearchCursor(r"C:\Users\cchen\Documents\temp\test.gdb\linetest")

point = arcpy.Point()
array = arcpy.Array()

featureList = []
cursor = arcpy.InsertCursor(r"C:\Users\cchen\Documents\temp\test.gdb\newlines")
feat = cursor.newRow()

for in_row in in_rows:
    # Set X and Y for start and end points
    point.X = in_row.sx
    point.Y = in_row.sy
    array.add(point)
    point.X = in_row.ex
    point.Y = in_row.ey
    array.add(point)
    # Create a Polyline object based on the array of points
    polyline = arcpy.Polyline(array)
    # Clear the array for future use
    array.removeAll()
    # Append to the list of Polyline objects
    featureList.append(polyline)
    # Insert the feature
    feat.shape = polyline
    cursor.insertRow(feat)
del feat
del cursor

It gives me this error message:
Runtime error
Traceback (most recent call last):
  File "", line 14, in
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\arcobjects_base.py", line 89, in _set
    return setattr(self._arc_object, attr_name, cval(val))
RuntimeError: Point: Input value is not numeric

Since the field is numeric(double), I am confused on why this does not work.


texts: Quantifying level of urbanization of catchment over time?

texts:
I am studying the impact of urbanization on the discharge levels of a river during flooding events. For this, I need to determine the catchment's percent urbanized area through land classification from the years 1992 to 2019.

I've thought of using Landsat images, fmask, ArcMap, eCognition, and the following procedure:

Use the Fmask algorithm for cloud cover removal.
Fill the removed portions using another (but temporally close) Landsat image with ArcMap.
Clip the final image with the catchment shapefile.
Have eCognition determine the percent urbanized area.

My problem is how to fill the portions that were removed by Fmask (step 2).

Is this procedure a viable method?

Moreover, I get confused as to how to process the Landsat images in ArcMap what with it having multiple bands.

And although I haven't really encountered any significant problem with eCognition aside from the somewhat low resolution image I get when working with my Landsat 8 data, I seek any suggestion as to how I can effectively use the program for determining the percent urbanized area.

I am relatively new at GIS.


texts:
Shape files of Europe maps are available on the Eurostat website. I downloaded NUTS_2013_60M_SH. The tool (readOGR from the R package Rgdal) I use to load the map shows the following proj4string:

 "+proj=longlat +ellps=GRS80 +no_defs"

Metadata provided with these files doesn't contain an EPSG code for the projection. I would like to use such a code to define the projection when I load the map.

A search for RS80 on the spatial reference site returns several code pages and none seem to be relevant.

What is the EPSG code of the projection used in these Eurostat maps?


texts: Use cell stats to loop through daily min/max rasters

texts:
I have a '2015' folder within which I have: 1) a 'Daily Min' 2) a 'Daily Max' and 3) a 'RasterCon' folders. 'Daily Min' and 'Daily Max' have shapefiles in them that I convert to raster and run a Con statement on and output to the 'RasterCon' folder. The Con rasters have values of 0 - 1.

In 'RasterCon', there are 2 files for each day (Day1ConMin, Day1ConMax, Day2ConMin, Day2ConMax). I want to run Cell Statistics sum function on each day's files - Day 1 Min/Day 1 Max, Day 2 Min/Day 2 Max. Ideally, the values should be 0 - 2. The last piece of my code is supposed to do this but does not return any output. Where am I going wrong?

import arcpy
import os
from arcpy import env
from arcpy.sa import *

#Set workspace to Max Temps folder and allow for overwrite
env.overwriteOutput = True
env.workspace = r"C:\Users\uma.bhandaram\Desktop\2015\DailyMax2015"
outputFolder = r"C:\Users\uma.bhandaram\Desktop\2015\RasterCon2015"

#Get Max Temp shapefiles and convert them to Raster
shpFileList = arcpy.ListFiles("*.shp")
for shpFile in shpFileList:
    shpFileName = os.path.splitext (shpFile)[0]
    rastFile = shpFileName + ".tif"
    arcpy.PointToRaster_conversion(shpFile, "Max_2015", rastFile, "MOST_FREQUENT", "", 0.05)

#Get raster files and run through Raster Calculator to see if Max threshold has been met
rastFileList = arcpy.ListFiles("*.tif")
for rastFile in rastFileList:
    rastFileName = os.path.splitext (rastFile) [0]
    MaxRstIn = Raster(rastFile)
    MaxRstOut = Con(MaxRstIn, 1, 0, "VALUE <=85")
    output_raster_path = os.path.join(outputFolder, "Con" + rastFile)
    MaxRstOut.save(output_raster_path)

#Set workspace to Min Temps folder and allow overwrite
env.overwriteOutput = True
env.workspace = r"C:\Users\uma.bhandaram\Desktop\2015\DailyMin2015"
outputFolder = r"C:\Users\uma.bhandaram\Desktop\2015\RasterCon2015"

#Loop through Min Temp shapefiles and convert them to Raster
shpFileList = arcpy.ListFiles("*.shp")
for shpFile in shpFileList:
    shpFileName = os.path.splitext (shpFile)[0]
    rastFile = shpFileName + ".tif"
    arcpy.PointToRaster_conversion(shpFile, "Min_2015", rastFile, "MOST_FREQUENT", "", 0.05)

#Loop through raster files and run through Raster Calculator to see if Min threshold has been met
rastFileList = arcpy.ListFiles("*.tif")
for rastFile in rastFileList:
    rastFileName = os.path.splitext (rastFile) [0]
    MinRstIn = Raster(rastFile)
    MinRstOut = Con(MinRstIn, 1, 0, "VALUE >=45")
    output_raster_path = os.path.join(outputFolder, "Con" + rastFile)
    MaxRstOut.save(output_raster_path)

#Loop through corresponding Min and Max rasters and calculate Cell Statistics (SUM)
#This will identify which days met both Min and Max thresholds (SUM = 2)
MainFolder = r"C:\Users\uma.bhandaram\Desktop\2015\RasterCon2015"
for (path, dirs, files) in os.walk(MainFolder):
for dir in dirs:
    env.workspace = os.path.join(path, dir)
    listRstrs = arcpy.ListRasters("Con.tif")
    for raster in listRstrs:
        rasterlist = []
        rasterlist.append(raster)
        outCellStats = CellStatistics(rasterlist, "SUM", "DATA")
        outCellStats.save(os.path.join(path, dir) + "Check.tif")


texts:
I new man in processing and modeling rasters by GRASS GIS. But I have one unsolved problem.

I have a lot of rasters with similar names: raster_1, raster_2, ..., raster_100.

Also in my work environment there are more rasters with unique names.

I need to do some steps:

seting to all rasters with names "raster_..." null values against zero;
calculation total (for all) minimum and maximum values for rasters with names "raster_...";
setting total for all rasters (with names "raster_...") color
palette based on total minimum and maximum values.

How can I to solve this problem? Do I must to use batch processing with names mask? Or GRASS give another ways?


texts: Script to extract UTM Zone using Range and Meridian information

texts:
I am looking for a python or php script that will accept the Range and Meridian values of a location and use it to identify the UTM zone.

I have a table that comprises of Easting, Northing, Range and Meridian but I'm missing the UTM zones of these points. I am able to manually identify the UTM zone knowing the range and Meridian of the location. I wish to automate this by using a python script that will scan in these parameters as variables and extract the UTM zones based on the parameters.

Can someone please help me out! It will really be exciting to achieve this as I really need to get it going, its so clumsy having to go through thousands of point data to this manually.

Meridian      Range      UTM Zone

   W1M       1 to  30    14

   W2M       1 to  30    13

   W3M       1 to  15    13

             16 to 30    12

   W4M       1  to 30    12

   W5M       1  to 30    11

Every data collected is described by a combination of meridian and range value. e.g NW-05-046-06-W4M to NW-05-046-05-W5M.


texts: Creating a Toolbox in ArcGIS Pro that can export results into an image format

texts:
How can I create a Toolbox in ArcGIS Pro that can export results from point and line feature classes and combine them into just one image format (e.g. .jpg, .png).
like so


texts: Choosing tool to add field based on nearest neighbor using ArcMap

texts:
Using ArcMap, I have two feature classes. The first is my polygon parcel layer with crop type, irrigation type, owner, etc. The second is my line canal layer, which has a field for our water schedule (how many days until a farmer can take water again).
I need to determine if our water schedule has an impact on what methods farmers use to irrigate their fields.
Using the ESRI basic license, what tool would I use to add the water schedule field from the canal layer to my parcel layer, based on the closest canal? The canals do not intersect the parcel boundaries.


texts: Error using two raster with different number of bands for Export Training Data For Deep Learning in ArcGIS Pro

texts:
I want to create a deep learning training dataset for pixel-classification which combines 3-band RGB data and 1-band elevation data. I use the ArcGIS Pro function "Export Training Data For Deep Learning". The input consists of the folder containing both rasters as well as the shapefile containing classified data. I'm able to create training datasets of multiple 3-band rasters, multiple 1-band rasters, but when I try to mix 3-band and 1-band I get following error:
"ERROR 002860: Tool parameters are inconsistent with the data you are trying to append to.
Failed to open raster dataset [Failed to open raster dataset].
Failed to execute (Export Training Data For Deep Learning)."
I assume the problem is either with the different number of bands, pixel depth or something else. Raster information for both datasets is shown below.

Here are the parameters for my function execution:

Where lies the problem and what would be the solution?


texts: GDAL import error, broken both python script and Qgis

texts:
I have both a standalone script trying to use gdal and Qgis both of which have run into an error with importing gdal:
My script errors:
line 18, in <module>
    import gdal
  File "/usr/lib/python3/dist-packages/gdal.py", line 2, in <module>
    from osgeo.gdal import deprecation_warn
ModuleNotFoundError: No module named 'osgeo.gdal'

and if I try from osgeo import gdal:
line 18, in <module>
    from osgeo import gdal
ImportError: cannot import name 'gdal' from 'osgeo' (/home/justinlarking/.local/lib/python3.8/site-packages/osgeo/__init__.py)

This is also affecting my Qgis as in I cant open it:

Full Qgis error:
Couldn't load plugin 'shapetools' due to an error when calling its classFactory() method

ImportError: cannot import name 'gdal' from 'osgeo' (/home/justinlarking/.local/lib/python3.8/site-packages/osgeo/__init__.py)
Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 334, in _startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/home/justinlarking/.local/share/QGIS/QGIS3/profiles/default/python/plugins/shapetools/__init__.py", line 10, in classFactory
    from .shapeTools import ShapeTools
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 792, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/justinlarking/.local/share/QGIS/QGIS3/profiles/default/python/plugins/shapetools/shapeTools.py", line 5, in
    import processing
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 792, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/__init__.py", line 35, in
    from processing.tools.raster import *  # NOQA
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 792, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/tools/raster.py", line 26, in
    from osgeo import gdal
ImportError: cannot import name 'gdal' from 'osgeo' (/home/justinlarking/.local/lib/python3.8/site-packages/osgeo/__init__.py)

Python version: 3.8.5 (default, Jul 28 2020, 12:59:40) [GCC 9.3.0]
QGIS version: 3.16.2-Hannover Hannover, f1660f9da5

Python Path:
/usr/share/qgis/python
/home/justinlarking/.local/share/QGIS/QGIS3/profiles/default/python
/home/justinlarking/.local/share/QGIS/QGIS3/profiles/default/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python38.zip
/usr/lib/python3.8
/usr/lib/python3.8/lib-dynload
/home/justinlarking/.local/lib/python3.8/site-packages
/usr/local/lib/python3.8/dist-packages
/usr/lib/python3/dist-packages
/home/justinlarking/.local/share/QGIS/QGIS3/profiles/default/python
/home/justinlarking/.local/share/QGIS/QGIS3/profiles/default/python/plugins/shapetools/ext-libs

Can anyone help me understand why my GDAL has broken? I did remove it and reinstall it however it is the same version as it was previously.


texts: Return rows after a spatial join where the count is zero

texts:
I am trying to determine the number points in a table of polygons AND include a row of zero for those instances where a polygon has no points in it. This query works well except that no zero rows are outputted:

SELECT COUNT(h.points) AS "num_points", ws.polygons
                 FROM points_table as h
                 JOIN polygon_table as ws
                 ON ST_Intersects(ws.geom, h.geometry)
                 GROUP BY ws.grouping_variable;

I tried using COALESCE like so:

SELECT COALESCE (COUNT(h.points), 0) AS "num_points", ws.polygons
                 FROM points_table as h
                 JOIN polygon_table as ws
                 ON ST_Intersects(ws.geom, h.geometry)
                 GROUP BY ws.grouping_variable;

This successfully ran but didn't return any instances where the row equal zero. I think I am getting a little confused about the order of SQL. I am quite certain that one route is to take the "non-zero" output and LEFT JOIN it with the original polygon_table but I am not clear where that should actually happen.

In summary - how do I return rows after a spatial join where the count is zero as well as all the other counts?


texts:
I have a DEM layer from WCS.
How can I download the data to process it to get slope, aspects, etc in QGIS?


texts: r.grow.distance with r.mask returns errors on a global scale, but not with a smaller extent or without a mask

texts:
We are trying to calculate distances to settlements on a global scale (with a 10 arc second resolution), to eliminate 2/3 of the area (water/ocean) to calculate we used land as a mask (r.mask) but this isn't working on this scale.

Without a mask it works fine but it takes 8 hours to run

With a mask we get the following errors:

WARNING: ZSTD compression error -14: Unsupported frame parameter
ERROR: Error uncompressing raster data for row 33288 of
CalledModuleError: Module run None ['r.grow.distance', '--o', '--q', '-m', 'distance=r101535547', 'metric=geodesic', 'input=r101535531'] ended with error

We checked the following:

Row 33288, there seems to be nothing strange in that row (there is data and all settlement are on land)
We used smaller datasets (with row 33288) with about 1/4 of the globe and with the mask and they all worked fine
Running it in ArcGIS with a mask, this worked fine

It runs on a E64s_V3 Azure machine with 432 GiB RAM reading/writing to a shared storage

Windows Server 2016 Datacenter

Grass 7.6

ZSTD 1.3.2.0 (the name of the dll in the OSGEO4W folder)


texts:
After I smoothed a layer of lines I lost the connection with some vertices.

Is there an automatic way to connect the various lines to avoid doing it manually?

I think the best thing is to use v.clean but I am not sure of the right parameters that should be set.