texts: Edit and query buttons disabled in geoexplorer

texts:
I have an issue: edit and query buttons are disabled in geoexplorer. This app use to be OK and then went wrong without doing anything on my server. How can I fix this ?

Thanks


texts: import arcgisscripting gives ImportError: DLL load failed: %1 is not a valid Win32 application

texts:
i am trying to import arcgisscripting and getting following error.

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
    import arcgisscripting
ImportError: DLL load failed: %1 is not a valid Win32 application.

I have windows 7 64 bits, i have installed python-2.7.3.amd64, pywin32-219.win-amd64-py2.7
and using arcgis9.3. i have set pythonpath accordingly which is set to python and arcgis/bin folder as well. But still getting the same error.

I have spent almost whole day to fix this issue and tried all the exisiting solutions but didnt get any success, i think my arcgis 9.3 is not 64 bits.


texts: rasters as independent variables in logistic regression: must they have the same values classification?

texts:
I am planning to build a predictive model for land productivity using arcGIS 10.1. and logistic regression. In my model, the Independent Variables would be, among others, distance from rivers, distance from spring-line, distance from the coastline.

I am wondering if it is sounder/useful/essential to use the same value classification for all the three aforementioned rasters, or if it would not pose any problem in analytical perspective to have rasters with different classification.


texts: Trying to re-project a folder of shapefiles using Python Script

texts:
I am trying to reproject multiple shapefiles in a folder. When I run my program this is the error I am getting:

line 8221, in Project
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000670: output Output Dataset or Feature Class is same as input Input Dataset or Feature Class
Failed to execute (Project).

I checked my code and I'm not seeing what Im doing wrong. What could be the issue here?

# Import system modules
import arcpy
import os

# Set environment settings
arcpy.env.workspace = "C:\\users\\data\\shapefile_roads"
arcpy.env.overwriteOutput = True
outWorkspace = "C:\\users\\data\\shapefile_roads"

for infc in arcpy.ListFeatureClasses():

        # Determine if the input has a defined coordinate system, can't project it if it does not
    dsc = arcpy.Describe(infc)

    if dsc.spatialReference.Name == "Unknown":
        print ('skipped this fc due to undefined coordinate system: ' + infc)
    else:
        # Determine the new output feature class path and name
        outfc = os.path.join(outWorkspace, infc)

        # Set output coordinate system
        outCS = arcpy.SpatialReference('WGS 1984')

        # run project tool
        arcpy.Project_management(infc, outfc, outCS)

        # check messages
        print(arcpy.GetMessages())


texts: Seeking tool to Convert Shapefile (.shp) to GTFS?

texts:
I need to convert a Shapefile (.shp), shx or dbf to  GTFS (General Transit Feed Specification https://developers.google.com/transit/gtfs/reference).

Do you know any tool for this issue?


texts: GDAL grid Not recognizing csv as having geometry

texts:
NON-DUPLICATE EDIT: This question is not a duplicate because I have tried the solution presented as the answer to the "original" question and it did not work. See results below.

I am attempting to use gdal.Grid in python to interpolate a .txt or .csv to a raster. I keep getting "ERROR 1: Failed to open datasource 'file.csv'. No point geometry found on layer file, skipping."

How do I get gdal.Grid to recognize the file? Or how do I need to alter the file in order to get gdal to recognize it?

EDIT: Ran ogrinfo twice, both with txt and csv versions. The results were the same each time, with the only difference being the extensions each time matched what was in the vrt and what the file was. ___Stupid question, if a text file is structured like a csv you can usually just change the extension to csv and programs will handle it, I've been doing that when "changing format", could this be part of the issue?

Below, column 1 is x, column 2 is y, column 3 is z but they do not have headers.

Data file structure:

555000,66600,001
555003,66610,002
555002,66601,003

VRT file

<OGRVRTDataSource>
    <OGRVRTLayer name="file">
        <SrcDataSource>CSV:file.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

Relavant python code

rasname = str(file)
rasname = rasname.split('\\')
rasname = rasname[-1]
rasname = rasname.strip('.txt')
vrtname = rasname + '.vrt'
vrtname = input_dir + '\\' + vrtname

outputfile = output_dir + str(rasname) + '.tif'

gridopt = gdal.GridOptions(format='GTiff',algorithm='invdist:power=2')

output = gdal.Grid(outputfile, vrtname, options=gridopt)

Results of ogrinfo

INFO: Open of 'file.vrt' using driver 'OGR_VRT' successful.
Warning 1: Cannot open datasource 'file.csv' in update mode. Trying again in read-only mode
Error 1: Failed to open datasource 'file.csv'

Layer name: file
Geometry: None
Feature Coutn: 0
Layer SRS WKT:
(unknown)


texts: ogr can't read shapefiles or json files even though geojson and esri shapefile are included in driver list

texts:
I am running gdal/ogr on mac os 10.12.4. I recently reinstalled gdal and now am having an issue when I try to read shapefiles or geojson files - I get the error: Unable to open datasource '/path/to/shp' with the following drivers. And in the driver list I can see ESRI Shapefile as well as GeoJSON. I tested reading a kml and it did not have issue with that. Does anyone know what might be wrong?

I also just checked my installations and it appears I have 3 installs of gdal - I know I have v1.11 and v2.1, not sure what the 3rd one but maybe this is causing issues?


texts: Import raster using raster2pgsql with out-of-db option not working

texts:
I am using raster2pgsql to import DTED data into postgis database.
I would like the data to be out-of-db, I am using the command as follows.

raster2pgsql -R /home/gis/elevation/DTED/area.dt2 public.dted

But the sql generated from this shows that the data is inserted into the table and there is no reference to the filename.
What am I missing here?

Please help.


texts: How to export raster into ascii without losing the statistics?

texts:
I had a problem in Mike Urban with a difference between a DEM (Digital Elevation Model) and a 2D overflood model and the support group requested me to remove the difference because the problems were at the edges of the area.

Then I found an inspiration by: ArcGis--how to lower a DEM's Z value along a line

I made a polyline in Arcgis and afterwards converted it to raster and reclassify the file(to remove the NoData) and at last I subtracted the area from the original DEM to edit the height by the Minus Tool.

Afterwards I tried to convert the raster file to ASCII by the converting tool, but when I open the file in MIKE Urban (GIS-based urban modeling system for water distribution systems and wastewater collection systems) the statistics (max, mean etc.) were missing.

Is there another way to convert the file without losing data?


texts:
I try to match two layers. First one is point-type with cities in Europe. Second one is grid layer (polygon) with elevation value. What I want to achieve is to to join the "Z" value from grid layer to point layer.
Both have the same CRS and the result is very strange.

I was trying to set different CRS for these layers but the result is still the same.

Only two "Z" values occurs in created layer.
Probably something is wrong with elevation data but I don't have any idea to solve it.
Here I attached these two layers
http://www.speedyshare.com/GVUqM/DEM.zip


texts: Why does "Join attributes by location" fail to compute sums when used in a model?

texts:
If I use "Join attributes by location" from the menu it works fine but if I'm using it within the Processing modeler it's not working.

I want to sum up previously selected records having a field area (calculated by "Export/add geometry columns").

Surprisingly, the sum fields are created but all records contain NULL.

Did I miss something? Any well known bug which I haven't found with Google? Do these two algorithm works differently?

BTW: If I run "Join attributes by location" from the processing toolbox, the algorithm works fine too.

(I'm using QGIS 2.8.1-Wien.)


texts: Load ASCII rasters in QGIS as singleband pseudocolor by default

texts:
I'm loading over 400 ASCII rasters into QGIS.
They all get displayed as 'singleband grey' by default.
Is there a way to change the style to 'singleband pseudocolor' by default?
I also want to set the min-max values in the TOC to the actual min-max values in the raster (all rasters have different min-max values).

I know you can do this layer by layer manually using

Layer properties
Style
Band rendering/render type: Singleband pseudocolor
Load min/max values - accuracy: Actual (slower)

But that's undoable for 400+ rasters.
So how can I set loading rasters to 'singleband pseudocolor' with actual min-max values by default?


texts: Reprojecting shapefile using DotSpatial gives unexpected lat/long coordinates?

texts:
I need to take an uploaded polygon shapefile and extract some information from it, including the lat/long positions of the centroids of the features.

Using DotSpatial I attempt to Reproject the shapefile into WGS 1984 to get latitude and longitude

if (!shapefile.Projection.IsLatLon)
{
    shapefile.Reproject(
        KnownCoordinateSystems.Geographic.World.WGS1984);
}

but the coordinates I get out land in South Dakota when they should be in Colorado.

Some centroids prior to reprojection:

{(3118274.90587493, 1458696.3665899)}
{(3118274.90587493, 1458696.3665899)}

The coordinates in South Dakota that I get after running the code above:

{(-104.040689162766, 43.4531513926421)}
{(-104.040690516737, 43.4529967532427)}

Am I missing a step?

EDIT:
Contents of .prj file
PROJCS["NAD_1983_HARN_StatePlane_Colorado_North_FIPS_0501_Feet",GEOGCS["GCS_North_American_1983_HARN",DATUM["D_North_American_1983_HARN",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",3000000.000316083],PARAMETER["False_Northing",999999.999996],PARAMETER["Central_Meridian",-105.5],PARAMETER["Standard_Parallel_1",39.71666666666667],PARAMETER["Standard_Parallel_2",40.78333333333333],PARAMETER["Latitude_Of_Origin",39.33333333333334],UNIT["Foot_US",0.3048006096012192]],VERTCS["NAD_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Meter",1.0]]


texts: How to get all Vector layers from OpenLayers map?

texts:
OpenLayers.Map object exposes a method call getLayersByClass.  How do I use this method to retrieve all of the Vector layers from the map?

I tried map.getLayersByClass("Vector") but this returned an empty array.


texts: Building extraction with LiDAR data (Improve Process)

texts:
I have the following process for automatically creating building outlines/footprints from Las files and was wondering if anyone can help improve it or has suggestions for other processes to use?

EXPERIMENTS:

A. ArcGIS - Lidardataset tools : while you can see just the building class points you can not create polygon automatically. There is the minimum bounding geometry but as it looks at the dataset as a whole it requires you to manually choose the points for each building.

B. External Lidar Analysis tool by Overwatch - building extraction is a core functionality. http://www.featureanalyst.com/lidar_analyst.htm (No longer active - July 28/2014)
It took approx 8 minutes to create a bare earth (essential) and then extract the buildings. It gives the best results by far, even capturing building under the trees (1&4). This is without even looking at a tutorial as it is a very intuitive interface.

Some issues are it has missed the covered corridors (3) but has correctly left out the open roof (2).

C. LASTOOLS - used lasmerge to merge 4 relevant tiles and then lasclip to extract zone of AOI only.

arcpy.lasboundary("AOI.las","1","buildings","true","true","shp","AOI_buildings.shp","","#","true")

The output was very jagged so I used:

arcpy.SimplifyBuilding_cartography("AOI","AOI_buildings_simplified.shp","4 Meters","0 SquareMeters","NO_CHECK")

and then ran:

arcpy.SimplifyPolygon_cartography("AOI_buildings_simplified","AOI_s.shp","POINT_REMOVE","1 Meters","0 SquareMeters","NO_CHECK","KEEP_COLLAPSED_POINTS")

1 and 2 show the issue areas. Part of it is due to bad lidar classification.

This gives good building lines for long straight edged buildings but if there is overlap by branches or it is a complex building with multiple edges, roofs on multiple levels close to each other etc, manual creation is necersary for any detailed site level work. For general estimation roof shape the process is accurate enough. Process takes less than 10mins total once the process has been finalized.

My Research: https://stackoverflow.com/questions/5037588/simplified-or-smooth-polygons-that-contain-the-original-detailed-polygon and the process suggested by @gromix
and What are Definition, Algorithms and Practical Solutions for Concave Hull?


texts: Joining attributes from nearest point in QGIS?

texts:
I have two layers with points, call them A and B, each with different attributes.

For each point in A I want to join all the attributes from the nearest point in B.

Ideally this would also include a field measuring the distance between both points.

I know how to so this in ArcGIS: right click layer A>joins and relates>joins, then join data based on spatial location and choose to join attributes from nearest point.

I could not find similar functionality in QGIS.  Any pointers?


texts: Viewshed analysis and coordinate systems ArcGIS 10.1

texts:
I am trying to run a viewshed analysis and I've downloaded DEM data from the UK ordnance survey (Landform-panorama). Since the data is in ASCII format I have created a seamless raster mosaic (and used the projected coordinate system - British National Grid). Now I am trying to display a point layer, with lat/lon coordinates of about 10 sites on top of the DEM raster layer. I have been told that both the raster layer and the point layer needs to have the same coordinate system, but when I use the projected coordinate system on the point layer, my points do not show up. However, if I use the geographic coordinate system on the point layer then the site locations pop up:

I was wondering whether any one could let me know if this is alright for a viewshed analysis? if not whether anyone could explain in layman's terms what I need to do in order to fix this issue. As always any advice would be greatly appreciated!


texts:
I'm learning Python and this has been a bit of a challenge!

My DataFrame looks like this below.

             name                   geometry
0     Uimastadion  POINT (24.93045 60.18882)
1  Hartwall Arena  POINT (24.92918 60.20570)
2           Talli  POINT (24.92607 60.21346)

I would like to create some sort of function on Python that takes an input (an address or name of a single point) then returns the nearest (say 10 0r 100) neighbors including their distances.

I've found this below but its not what I'm looking for.

https://automating-gis-processes.github.io/site/notebooks/L3/nearest-neighbor-faster.html

How can I do this?


texts: Determining distance on line of point, given measure?

texts:
I need to determine into a standalone python script the distance of a point on a polyline with M-values and its' coordinates, at a specific measure.

For example, I have to calculate the X,Y of a point situated at measure 12.345 on a polyline with the measures from 10 to 15 and considering that the startpoint of the polyline has m=10

What I did so far was to use the .positionAlongLine method, with the distance= 12.345-10, which in field calculator of the polyline featureclass would look like:

!Shape!.positionAlongLine(2.345,False).Centroid.x

The problem is that the resulted point with the calculated XY is not exactly at measure 12.345 and not at the distance 2.345 towards the startpoint, as expected.


texts:
The Question:

Is there way to leverage GIS tools to author/edit a schematic map where elements' geometries has canvas-based coordinates rather than real-world coordinates on Earth?

Perhaps to achieve this, a special Coordinate Reference System (CRS) could be used; one that does "No Projection", one that maps an (X,Y) point directly as-is on a Canvas (instead of on Earth).

Motivation:

I'm maintaining a Geo-database feeding a real-world map. I'd like to also maintain, in parallel, a "schematic map" where I arbitrarily locate the same elements on a planner canvas. So while each element has a Geometry column with real-world coordinates, it also has another geometry column (point, line...etc) with simple X/Y canvas-based coordinates (not real-world) (similar to vectors in an SVG).

Then it would be of great value if I could view/edit this schematic map in GIS tools (like QGIS).

  Example:

  A Real-World Map, in one QGIS project, has layers of features using their typical geometry column in a typical CRS (say WGS 84) - this map is typical, no problems here!
  A Schematic Map, in another QGIS project, has layers of the same features but using another geometry column that is canvas-based (rather than the real-world) coordinates.

Although the essence of GIS is to manage elements located on Earth (real-world coordinates), my special use-case would benefit highly from managing a schematic map in the same environment and with the same tools as a real one.

Any thoughts?

NOTE: I'm not trying to generate a schematic map automatically (like ArcGIS' Schematics extension does).


texts: Model Builder - Integrating Select by Location and Buffer

texts:
So I'm trying to create a model tool in ArcGIS 10.1 that will integrate select by location and buffer.

These are the features that the user is supposed to be able to choose:

Select Layer by Location - selecting features

Buffer - distance of buffer radius (must be in km)

These are the features that the users will not be able to change:

Select Layer by Location - Relationship (COMPLETELY_WITHIN), Input Feature Layer

Buffer - Dissolve type (LIST), End type (ROUND), Side type (FULL)

This is the model I have created so far, but it doesn't let me enter all the fields I would like to enter, and when I run it it generates no output:

How do I create a model where I can enter the values I'd like?


texts: Draw rectangle with a given size (height/length) from attribute table

texts:
I would like to build a rectangle (shape or symbol) with a given size (height/length) from the attribute-table.

In detail:
I build a Atlas of maps with alternating sizes, like 145x129 or 165x129. The size of every feature is written in the attribute table of the coverage-layer. So far so good: every map got his own defined size using data defined override for map length and height.

But now I would like the bounding boxes of the atlas-maps to be shown as overlays on the map in my canvas. (Like you can do with the overlay-function in the print composer)

The best approach seemed to be to use the geometry-generator and build a polygon around the centroid of the feature. I already started a post months ago: Using QGIS Geometry Generator to get rectangle from point?
With the result, that it is not possible because the "Geometry generator does not support calculations inside WKT".
I still have that problem to solve. Do anyone know a other approach?


texts:
Have been playing around with this answer for a while and I can't get it to work in QGIS 3.4.4:

https://docs.mapbox.com/help/tutorials/mapbox-arcgis-qgis/#add-mapbox-maps-in-arcmap

No help on YouTube or MapBox Tutorials for this scenario.

MapBox keeps saying that works for QGIS 2.x - si that something to do with it?

Tried import via WMS/WMTS (Doesn't connect), via XYZ tiles (raster link shows connection is fine, but doesn't display), via OWS (Doesn't connect).


texts: Using ArcGIS Field calculator with global variable?

texts:
I have an issue with creating a field calculator expression in python that utilizes global variable.

How my table looks like:

This is a fragment of road network shapefile. Each road has a code (numer), but some segments of the network have more than one code. I need to duplicate the segments that have multiple codes (already done, OrygFID is the FID of polyline from before duplicating it), assign one code per one polyline and dissolve them by their code - so I get one polyline per road.

This is how my field calculator expression looks like:

count = 0
globfid = ""

def function (numer, fid):
 global count
 global globfid

 codes = numer[3:-1]
 codes_s = codes.split(",")

 if fid == globfid:
  globfid = fid
  count += 1
  return codes_s[count]
 else:
  globfid = fid
  return codes_s[count]

It's supposed to clear numer from everything but the road codes (works), separate the codes into a list (works). Then it's supposed to check for OrygFID of previous row (stored as global variable), and if it's the same then it should return next code from the stored list (doesn't work). I call it using function (numer, OrygFID). The result I get for the above table is as follows (one per row).

NumerDrg: 1153K, 1147K, 1155K, 1155K.

How it should look like: 1153K, 1147K, 1155K, 1156K.

I am extremely new to Python and kind of can't wrap my head around this. Any idea why my script won't work?


texts: Removing map from ArcGIS Pro project using ArcPy?

texts:
I have a desirable requirement to be able to remove a map from an ArcGIS Pro project using ArcPy.

However the ArcGISProject class does not seem to provide a method to do this.

Is it possible to do so by any other means from ArcPy?


texts:
I am trying to copy feature layers using the following code, but get no results. Can anyone help me understand why the code does not work?

import os
import arcpy
arcpy.env.workspace = ws = 'C:\\WorkSpace\\sahie_saipe\\sahie_saipe.gdb'
dictionary_Layer = {'co06_Layer':'co06', 'co07_Layer':'co07', 'cnty08_Layer':'co08', 'co09_Layer':'co09', 'co10_Layer':'co10', 'co11_Layer':'co11', 'co12_Layer':'co12', 'st_Select':'state'}
for layer in arcpy.mapping.ListLayers(mxd):
    if layer in dictionary_Layer:
        old = os.path.join(ws, layer)
        new = os.path.join(ws, dictionary_Layer[layer])
        try:
            arcpy.CopyFeatures_management(old,new,"")
            print 'copy %s to %s' %(layer, dictionary_Layer[layer])
        except:
            print 'Could not copy %s'


texts: Converting standard to feature linked annotation?

texts: Is there a way in ArcGIS 10.2.2 & ArcSDE 9.3 to convert standard annotation to feature linked annotation?


texts:
I am trying to trace in QGIS after using GEO REFERENCER on the Vector menu of QGIS  to import an image . QGIS  which is amazing the youtube tutorials
 are brilliant videos are sadly for the 2.8 version and not 3.10.1

Using the earlier version videos  I have imported the image I have and its amazing , however when creating a feature on it as a polygon using 'Add polygon ' tool box icon. I can draw around what I want but only exit this polygon creation  by clicking 'pencil' tool box item . This creates a red  outline on the map . Now if I try and edit this like my previous GEOJSON layers  I cannot and it seems locked . Am I doing something wrong as I would rather just use QGIS than create in another product.

To explain :
I have created a standard 'New Shapefile' layer and have selected this and drawn a shape into the layer using 'Add polygon' tool. It appears as a red outline  which looks great but does seem to be a feature.IF i ask the layer to count the feature it shows (0) features  The second video on youtube for QGIS 2.8 shows you how to do this addition in 2.8 not 3.10.1 . Once this is in, I expected to be able to use the 'Vertex' Tool to edit points, however when I clicked on the outline it appeared locked and did not  change to a polygon with points . I am unsure if I have missed a step as I tried to save this to GeoJson text file and the file was empty . It would possibly appear then that this is not saved into the layer as the geojson files are of 0 bytes? I must not be doing something correctly as this is an awesome package .

There must be a difference between these GEOJSON files and the layers I am now creating?


texts: Map canvas saveAsImage() method saves wrong image in QGIS

texts:
I've been trying to write a small python script to output 2 tif's  based on given box extents.
The script all works, box coordinates convert to project coordinates, zooms to right place (would like it to pause between zoom changes), just the 2 outputs are wrong, they both show the image that was displaying before the script is run.

import qgis.utils
import time

source_crs = QgsCoordinateReferenceSystem(4326)#heographic
dest_crs = QgsCoordinateReferenceSystem(3857)
transform = QgsCoordinateTransform(source_crs, dest_crs)

myBox = QgsRectangle(0, 51,2,53)
new_box = transform.transformBoundingBox(myBox)

canvas = qgis.utils.iface.mapCanvas()
canvas.setExtent(new_box)
#canvas.clearExtentHistory()
#canvas.clearCache();
canvas.refresh()
canvas.saveAsImage("c:/data/qgis_image.tif" ,None ,"TIF")

myBox = QgsRectangle(-4, 51,0,53)
new_box = transform.transformBoundingBox(myBox)

canvas = qgis.utils.iface.mapCanvas()
canvas.setExtent(new_box)
#canvas.clearExtentHistory()
#canvas.clearCache();
canvas.refresh()
canvas.saveAsImage("c:/data/qgis_image2.tif" ,None ,"TIF")


texts: Selecting specific rows and column from a shapefile using Fiona and Shapely

texts:
I am new to using PyCharm, but I am trying to

cut out 4 rows from a county shapefile
from that 4 rows specify the columns or attributes to also be extracted
and output the result as a new shapefile.
but when I try this code, it gave me some errors, please any ideas?

SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 2-3: truncated \UXXXXXXXX escape

new error after rectifying file path
SyntaxError: positional argument follows keyword argument

import gdal_workaround
import fiona
import shapely

if __name__ == "__main__":

with fiona.open('C:\\Users\\Desktop\\GEM\\20test\\counties\\counties.shp') as input:

  output_schema = input.schema.copy()

  with fiona.open(C:\\Users\\Desktop\\GEM\20test\\counties\\counties2.shp , 'w', 'ESRI Shapefile',output_schema, crs=input.epsg, 4326) as output:
    for elem in selection:
        output.write({'properties': elem['nut2name'],'geometry': elem['polygon']})


texts:
I'm a bit new to ESRI and GIS so forgive me if this is a simple question.

Looking at this example ESRI Javascript sample (source code is here) I see you can click on a shape (via a published shape file) and if the point upon which you click is within a shape an action fires.

Is it possible to send a lat/long via some mechanism and IF that lat/long is within the shape the same (or similar) action fires?

or

Can I see if a lat/long is contained within ANY polygon that is contained in a published shapefile.

If my understanding of the ESRI javascript API is correct the point that the mouse is clicked does translate to a lat/long location.


texts: Is this Intersect or Union algorithm associative?

texts:
Similar to Is the Intersection Operation Commutative? . Spoiler for that one: Intersection is commutative and associative.

I've created an algorithm (in PostGIS but the theory would work in other languages) that takes an array of polygons as an argument and performs the following: if the polygons overlap, take their intersection. If they do not overlap, union them instead.

     i integer;
     tmpGeom geometry;
begin
     tmpGeom := geoms[1];
     FOR i IN 1..array_length(geoms,1) LOOP
        IF ST_CROSSES(tmpGeom,geoms[i]) THEN
              tmpGeom:= ST_Intersection(tmpGeom,geoms[i]);
        ELSE
              tmpGeom:= ST_UNION(tmpGeom,geoms[i]);
        END IF;
     END LOOP;
return tmpGeom;
end;

My question is, does the order of the operation matter? I can test that by randomizing the order for an input of a set of polygons, but I'm looking for some resources to prove it.


texts: Preventing flow accumulation in one watershed from traveling up tributaries into adjacent watersheds?

texts:
I am using TauDEM in ArcMap 10.3.1 to try and calculate flow accumulation and generate stream networks for HUC10 watersheds. I am using 3m DEMs acquired from the National Map.

The HUC10 I am working with consists primarily of a segment of a large river, and there are several large tributaries flowing into this river segment, but each tributary is part of a separate HUC10. You can see this in the following image, the HUC10 I am working with is highlighted in blue and the tributaries are coming in from the top of the picture:

I am trying to calculate flow accumulation only for the highligted HUC10, so I have clipped it from the rest of the DEM for my analysis. My specific problem occurs when calculating the flow accumulation (based on D8 flow direction) - the flow accumulation will start at the furthest upstream portion of the main river (left hand side of the above picture), and the raster cells will increase in value until they reach one of these incoming tributaries. At this point they intersect the edge of my clipped DEM as if it were an outlet point instead of continuing down the main river.

I have tried buffering the DEM various distances outward from the watershed boundary to try and avoid this problem where the flow accumulation intersects the edge and this works sometimes but other times it just allows the flow accumulation to flow further up the incoming tributaries. See photo below, the flow accumulation grid (based off of a DEM with a 500m buffer around the watershed boundary) is displayed at 50% transparency on top of an aerial photo for reference:

Does anybody have any ideas on how I can make the flow accumulation continue down the main river instead of up the incoming tributaries?

I'm sorry if my explanation of the problem was not adequate, I have never used the GIS stack exchange before. Please let me know if you need any further information or clarification.


texts: Problem with automated data storing during raster::crop (‘legend’ is not a slot in class “RasterStack”)

texts:
I am working with a set of raster layers containing CHIRPS monthly rainfall amounts over Africa. For purposes of reproducing the below code, I staged the decompressed tif files on
Dropbox.

Now, let's say I want to crop the stacked images by a given geographic extent.
Here's the referring code.

## file import
library(raster)

ch_fls_chirps <- list.files("path/to/tif", pattern = ".tif$", full.names = TRUE)
rst_chirps <- stack(ch_fls_chirps)

## desired extent (xmin, xmax, ymin, ymax)
ext <- extent(c(37, 37.7, -3.4, -2.84))

## cropping (with output storage)
rst_chirps <- crop(rst_chirps, ext, snap = "out",
                   filename = "path/to/tif/CRP",
                   bylayer = TRUE, suffix = names(rst_chirps),
                   format = "GTiff", overwrite = TRUE)

When performing raster::crop, the following error message is thrown:

Error in checkAtAssignment("RasterStack", "legend", ".RasterLegend") :
‘legend’ is not a slot in class “RasterStack”

By contrast, the code works just fine when output storage is omitted during raster::crop and manually done afterwards.

## cropping (without output storage)
rst_chirps <- crop(rst_chirps, ext, snap = "out")
rst_chirps <- writeRaster(rst_chirps, filename = "path/to/tif/CRP"),
                          bylayer = TRUE, suffix = names(rst_chirps),
                          format = "GTiff", overwrite = TRUE)

Can anybody tell me what the issue with direct output storage from within the
crop command is all about?


texts: How to cast geometry to varbinary in PostGIS / PostgreSQL

texts:
In MSSQL Server, casting geometry to varbinary(max) produces a string that looks like this:

0x3D0B0000010CA0F3BE3B5CE9474100EC414E3AF93941

When you insert that string into a geometry column, a geometry object is built.

However, I can't seem to find out what the equivalent data format is in PostGIS / PostgreSQL in order to bring that same varbinary format out of a PostGIS geometry column (for the purpose of loading it into an MSSQL Server Geometry object).

I've read that the 'bytea' format does the same as varbinary, but the bytea version of a PostGIS geometry object looks like this:

\001\001\000\000 \346\020\000\000\263\234\021o\224?Z\300-\265\373^\322\334C@

Is there a way to get the same string that the MSSQL varbinary format produces from PostGIS / PostgreSQL?

Some background into this: I'm developing my own ETL tools using Python.

psycopg2 works well with PostgreSQL (PostGIS) database conversion - and PostgreSQL casts geometry to varchar in a SQL statement in my python script, which inserts nicely back into a PostGIS geometry column - PostGIS ETL done!

pypyodbc is working great for the receiving end of a psycopg2 query, but as the problem above states, I can't find the correct format to cast the PostGIS geometry to in the SQL in the python script... otherwise, I can read other columns from PostgreSQL and write them to MSSQL just fine.


texts: Landsat8 View Imagery Service in ArcMap Layer Time Issue

texts:
I'm attempting to work with one of the newer ArcGIS Online subscription imagery services called Landsat 8 Views which allows you to view Landsat data for the US across several dates in 2013.  I can get the imagery service to click thru the various dates in arcgis.com viewer but not in ArcMap 10.2  When I attempt to enable time on the Landsat 8 Views layer in ArcMap it removes the layer from even being visible.  Any ideas what may be going on here and what may be wrong?  It's a huge amount of data in the attribute table and possibly part of the problem.  I would like to clip specific dates for small areas but I can only get a single date to show which I guess is the one displaying on top.

The server address for the Landsat 8 Views is here.


texts: Python and gdalwarp: change resolution of raster

texts:
I am new to GDAL and have been having some hard time figuring out how to change a raster's resolution and how to align it to another existing raster.

I have come to the point where I have a raster made of a numpy array class. I exported it into an ASCII file. It looks similar to this (as an example, this is the ascii file)...

test.txt:

ncols 4
nrows 4
xllcorner 20
yllcorner 8.5
cellsize 0.5
nodata_value -9999
0.1 0.2 0.3 0.4
0.2 0.3 0.4 0.5
0.3 0.4 0.5 0.6
0.4 0.5 0.6 0.7

As you can see, it contains no projection information.
On the other side I have an existing raster saved on the computer which is in WGS84 UTM32N, and has a different resolution.

I tried to google around to find how I can use GDAL to change resolution of the ASCII raster while aligning it to the existing raster on my computer.

Does the ASCII raster need to have some projection data or other kind of information? Could somebody explain to me how to use gdalwarp with Python?


texts:
I am working with ecological niche models, I have a raster which represents the current scenario (img). I have a raster  which represents the future scenario.
To work on R, analyzing the lost areas and the areas gained in the future I need to reclassify the rasters.
I need to traform the raster to binary, assigning a value of 0 to unsuitable areas, a value of 1 to suitable areas.
How can I do?


texts: Using Not Equal to Wildcard in ArcMap Definition Query?

texts:
I'm attempting to use a definition query in order to display only a subset of data.
I have a large CAD drawing that I have imported into my map. I do not need many of the lines from the .dwg file to be shown. I noticed that most of the CAD lines I did not want to display shared similarly named attributes. Some lines were "Poly-Line-1", some were "Poly-Line-3", "Poly-Line-A", so on and so forth.
The query I attempted to use was the following:

"Attribute" <> 'Poly-Line-%'

I was able to apply this statement, which seemed to indicate that it was syntactically correct. Unfortunately the file was not affected and each line was drawn as it had been before. I then tried:

"Attribute" <> 'Poly%'

Again, it was treated as if the syntax was correct, but no change took place. As a workaround I ended up using a query similar to

"Attribute" NOT IN ('Poly-Line-1', 'Poly-Line-2', 'Poly-Line-3', etc)

This was more efficient than something like

"Attribute" <> 'Poly-Line-1' AND "Attribute" <> 'Poly-Line-2' AND...

However it still took a while to list every thing I wanted.

I'm interested in any solutions that are more efficient than my final workaround and don't make use of ArcPy scripts. The problem I described is something I encounter at the workplace and I am mostly restricted to using the ArcMap GUI there.


texts: Leaflet markers from CSV with multiple columns employing omnivore

texts:
Goal: Show markers for lat-lon positions from a CSV file. The other columns need to be placed in pop-up as attributes.

CSV file can be seen at: airport.csv

Code:

headers = ['iata', 'title', 'city', 'state', 'country'] // determined programmatically and excludes latitude longitude

customLayerLocalData = L.geoJson(null, {
    onEachFeature: function (feature, layer) {
        var popupContent = '';
        $.each(headers, function (h) {
            popupContent += feature.properties + '.' + headers[h];
            if (h != headers.length - 1) {
                popupContent += '+'
            };

        });
        layer.bindPopup(popupContent);
    },
    pointToLayer: function (feature, latlng) {
        return L.circleMarker(latlng, customMarkerOptions);
    }
});
layerLocalData = omnivore.csv.parse(fr.result, null, customLayerLocalData).addTo(map);

Problem: The pop-up shows the object:

 [object Object].iata
+[object Object].title
+[object Object].city
+[object Object].state
+[object Object].country

If I change: layer.bindPopup(popup content); to layer.bindPopup(eval(popup content));, it works. But, I am not happy with this solution (does not feel right). Morever, in this case new lines can't be inserted inside the custom geoJSON layer.

I'm looking for a smarter solution that can interpret the 'popupContent' and pass onto to layer.bindPopup correctly.


texts: Using arcpy.da.UpdateCursor to updateRow gives RuntimeError: row contains bad value?

texts:
I have a fishnet grid (polygon SHP), and all the little fishnets have a "Zone_ID." There are generally several dozen to a few hundred fishnet polys associated with each Zone_ID. I'm trying to number each fishnet poly based on the zone. Screenshot below is a simple example of what I'm tring to code. (Fishnets have already been intersected with zones, and the topology is correct.)

arcpy.AddField_management("MyFishnet_wZones.shp", "Fishnet_ID", "SHORT", field_is_nullable = "NULLABLE")

nPolygons = int(arcpy.GetCount_management("MyZones.shp").getOutput(0)) # length nPolygons is the total number of zones

for n in range(0, nPolygons):
    counter = 1
    with arcpy.da.UpdateCursor("MyFishnet_wZones.shp", ["Zone_ID", "Fishnet_ID"]) as cursor:
        for row in cursor:
            while row[0] == n + 1: # because Zone_ID starts at 1, not 0
                row[1] = counter
                cursor.updateRow(row)
                counter += 1

And here's the error:

    File "", line X, in <module>
        cursor.updateRow(row)
RuntimeError: The row contains a bad value.


texts: Using ArcGIS Desktop Advanced to host map services?

texts:
I'm fairly new to ArcGIS but ben around GIS / mapping a long time. Part of getting up to speed will involve building map services and consuming them, etc, etc.

Can I use ArcGIS Advanced Desktop to host map services? This would be for learning only (at home so no access to ArcGIS server) so no need to be concerned with performance. At most, one other person might take a look at what I'm doing. If this can't be done, any suggestions for alternatives?

ArcGIS 10.1 now and 10.3 in the near future.


texts:
I executed a  query to know where a MultiLineString touches itself using st_touches from postgis as shown here :

select a.id,b.id ,a.geom,b.geom, ST_Touches(a.geom, b.geom)
from pipe AS a , pipe AS b
where b.id > a.id and ST_Touches(a.geom, b.geom) = true

yet the result is boolean (true or false).

I want to generate a point layer in each touching location. I tried to do create as select but it returned a table without geometry column.

How can I generate the point locations?


texts:
When i try to get an automatic fieldname with the following expression: coalesce("Borserie",NULL,'Skriv inn ledig nummerserie') it doesn't get the value/string from the table? The weird thing is that it does look right in the preview?

The table looks like this:

And the result looks like this:

Tried and tried for a couple of days but nothing seems to work out. Anyone with some good hints or a solution?


texts: Create a new grid layer from a CSV, given the four corner coordinates


texts: Applying gdalwarp cutline to PostGIS data

texts:
I can not figure out how I directly can read a polygon into gdawarp's -cutline argument.

My current approach is:

ogr2ogr -f "ESRI Shapefile" mydata.shp PG:"user=myuser dbname=gisdb" -sql "select * from polys where id='12345'"
and then I run

gdalwarp -cutline mydata.shp in.tif out.tif

How could I skip the step with the shapefile?


texts: Understanding *coords meaning when used with ArcPy?

texts:
What does *coords in the following code mean?

import arcpy

feature_info = [[[1, 2], [2, 4], [3, 7]],
                [[6, 8], [5, 7], [7, 2], [9, 5]]]

features = []

for feature in feature_info:
     Create a Polyline object based on the array of points
     Append to the list of Polyline objects
     features.append(
         arcpy.Polyline(
             arcpy.Array([arcpy.Point(*coords) for coords in feature])))

arcpy.CopyFeatures_management(features, "c:/geometry/polylines.shp")


texts: How can I symbolize my rasters in ArcGIS to match levelplot in R?

texts:
I need to create more customized maps in ArcGIS than what I have in R, however I can't seem to get the outputs to match. I have model outputs as ascii that I have converted to tiffs in R and plotted:

ras<-function(x,b){
  m<-x
  crs(m)<-"+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

  if (require(rgdal)) {
    writeRaster(m,filename=paste0(b,".tif"),driver="GTiff",overwrite=TRUE)
  }
}
ras(raster,"Example")
#plot
col<-rasterTheme(region=brewer.pal(8,"YlOrRd"))
levelplot(raster,par.settings=col)

Using those tiffs, I am trying to display the tiffs in ArcMap, however the surfaces look slightly different. Here is a screenshot to illustrate this:

On the right is the same raster in R using levelplot however there is a slight gradient outside the "hotspots" with values around 0.1-0.2. The same cannot be seen on the raster to the left.

Is there a way to get the symbology consistent across the two? I'm wondering if my conversion is the reason as to why these are different.


texts: Performing advanced labeling using Python in ArcMap?

texts:
I am attempting to write a Python label expression that would involve multiple label definitions.

Here is what I am trying to do using ArcMap 10.1.

I have a point shapefile representing Monitoring Wells. For each Well I have taken concentration readings for various chemicals. What I need to do is label all of these concentrations for each well and color code them based on whether or not they exceed a certain limit.

This is what I have:

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

How do I repeat this so I get a label that looks somthing like this:


texts: ArcPy: How best to handle output that could be a shapefile or feature class?

texts: In many Esri examples, including this one, they code in the output ".shp".  What is the best way if you want to be able to handle a shapefile OR a feature class output.  Let say if a GDB is provided as the workspace, then features class, if a folder directory then .shp.  Is there any ArcPy built-in to handing this?


texts: Adding xy data to table data from collector in ArcGIS Desktop?

texts: I can easily download all the data from the collector app to my PC, but how do you add the xy (coordinates) to be included as well?


texts:
I have been using the following python code on QGIS3 (2.99) to look up fields, for several months

id = feature.attributes()[layer.pendingFields().lookupField('ID')]

And it gives me an error today:

AttributeError: 'QgsVectorLayer' object has no attribute 'pendingFields'

Is this another breaking change in the QGIS development version? And what's the correct way to write it now?


texts: How to create a matrix of points inside a polygon with QGIS?

texts:
I need to create a matrix of points (e.g. 3.5 x 0.5 meters), with a defined orientation, inside a polygon. Essentially, I have to plan the arrangement of plants in a grove in order to plan an irrigation system.

I appreciate any hints

Antonino


texts: Is there an international standard for region codes that defines the European Union?

texts:
I have seen the ISO 3166 defines 2-letter, 3-letter and numerical codes for countries around the world. However, the European Union does not seem to have a code.
Is there an international standard for region codes that defines the European Union and maybe other political/economic unions/regions?


texts: Merging linestrings in QGIS using PostgreSQL and PostGIS

texts:
I've to merge all linestrings from a PostgreSQL table "river" where not more than two start-/endpoints are overlapping. In the illustration below I use a "M" to mark all linestrings that needs to be merged.

Unfortunately I'm new to PostGIS. To my knowledge QGIS needs an id (gid) and a geometry (geom) column to load a table or view. But I'm not sure which PostGIS function will do the job. Can anyone please help me with the SQL query?


texts: How to get feature id by using QgsExpression?

texts:
I want to get feature id of the attributes by using QgsExpression. I use this code:

layer = iface.activeLayer()
expr = QgsExpression('$id')
it = layer.getFeatures(QgsFeatureRequest(expr))
values = [i.id() for i in it]

However, in my layer, there are 20 points. It gives me 19 points. How can I fix the problem?


texts: Configure the geoserver gdal plugin on Debian buster running open Jre11

texts:
I am trying to setup the gdal plugin in a geoserver docker container for geoserver 2.16 series. I have installed libgdal java and gdal but I get the following error
No gdaljni in java library path.
  I have set the LD_Library path and downloaded the image io extensions and setup gdal_data env variables.
Most of the articles I have seen reference old links which are not useful.


texts: How to merge / join line segments by attribute and split at intersection

texts:
I have a complex network of line features with a lot of attributes (and unique IDs). Many features have the exact same values (exept the ID). I need to join all features with the same attribute which are next to each other.

that's what I have: an example network, many features with same values and attributes:
That's what I Want: the lines are merged together.

I want to merge them by topology: a line should be joined with a touching one. But: where more than 2 are touching (intersection) they should not be merged.

I tried:

dissolve, but then they are all joined, even when far away. then I tried to split up with "multi to single part" but then I came back to where I started.
Join Lines does what I want but not automatically and it's too much to select and join them by hand

Ideal would be a function that merges by 1.aligned neighbour (if more than 1, e.g. at intersection: no join) and 2. by same attribute value

It is very similiar to posted problems:
Join line segments sharing a common attribute, using GDAL
;
How to join segments of a line in QGIS?

But those have not been answered to a point where the issue is solved.


texts: Start and end point of line segments (line shapefile) and getting specific labels from a point shapefile

texts:
I have a line shapefile where I need to break the line at the specific intersection with points using QGIS 2.18.20. I refer to these segments as ‘tray sections’.

Now I need to have two fields: f_from and f_to in the ‘tray sections’. In these fields, I need to calculate the specific label from the point layer ie a specific ID number.  I also need to have x and y fields for the f_from and x and y for the f_to as well.

Eg. See screenshot below:


texts: Recommendations for a fre OSM based walking tour mobile app?

texts:
I am working with a small municipality with a very small (read:no) budget.  They would like to build a walking tour for mobile use based around natural features, and some history but mostly nature walks.

I have been looking around at the array of services and approaches to doing this (online serving of data vs. download in advance, media rich vs. just text, etc.) and so far have not found "the one".

I am hoping some folks here have had some experience building walking tours from easy-to-approach apps or services, as opposed to having an ArcGIS Online account, or building my own city-specific app which seems to be the way a lot of places have gone.

The lo-fi approach I want is something where I can upload a line feature for the route and point feature class for the stops, as KML or shapefile, and have that information overlaid on Open Street Maps for navigation.  The user would navigate to the points and tap the marker to read the text about the area.  Even better would be a location-aware app that allows multimedia and when the user walks within a proximity of the point, it notifies them that there is text/audio/images available to view.

Any recommendations for accomplishing this?  There seem to be a lot of options, so looking for simple workflow, and not really wanting to start downloading a myriad of apps to try them all out.


texts: Postgis Query to join based on two individual table to matching two individual table with value being seperated by 'comma" in both table

texts:
I have polygon(properties) in postgis table and have a info of excel which need to join on posting table property.  But it's not simple join it based on two column with matching individual column with separated value by "," .

Postgis table

Excel table

I want match if any of comma seperated match the value in respective column. In above case all should be matched.

So far I tried this script but not work good.

SELECT tbl1.gid, id,tbl1.matchid1,cast(tbl2."matchid1" AS varchar) ,  tbl1.matchid2,cast(tbl2.matchid2 AS varchar),  tbl1.geom
    FROM public.tbl1 tbl1 left join  tbl2 on (cast(tbl1.matchid1 AS varchar) ~~ cast(tbl2.matchid1 AS varchar)) or
    (cast(tbl1.matchid2 AS varchar) ~~ cast(tbl2.matchid2 AS varchar))


texts: Exclude geometry column from properties of ST_AsGeoJSON output

texts:
The ST_AsGeoJSON function with the following query:
SELECT ST_AsGeoJSON(t.*) FROM vw_ph t (ph, capture_datetime) LIMIT 1;

Returns the following:
  {
     "type":"Feature",
     "geometry":{
        "type":"Point",
        "coordinates":[
           -88.987035,
           30.161476
        ]
     },
     "properties":{
        "id":1,
        "capture_datetime":"2020-08-26T19:51:05.575-05:00",
        "ph":16,
        "boat_datetime":"2020-08-26T19:51:05.538-05:00",
        "latitude":30.161476,
        "longitude":-88.987035,
        "an_point_geom":{
           "type":"Point",
           "coordinates":[
              -88.987035,
              30.161476
           ]
        }
     }
  }

I am trying to get a geojson output without having the geometry column an_point_geom repeated in the properties field as below as it will add unnecessary size in a large query and polygon geometry.
The closest I have been to this is with the query below.
SELECT ST_AsGeoJSON((t.ph , t.capture_datetime, t.an_point_geom, t.boat_datetime)) FROM vw_ph t LIMIT 1;

However, this returns the result without the geometry column in properties but the columns are named f1, f2... which is not desirable.
  {
     "type":"Feature",
     "geometry":{
        "type":"Point",
        "coordinates":[
           -88.987035,
           30.161476
        ]
     },
     "properties":{
        "f1":16,
        "f2":"2020-08-26T19:51:05.575-05:00",
        "f4":"2020-08-26T19:51:05.538-05:00"
     }
  }

Is there any way to exclude the geometry column from properties without another query?
I am using PostGIS 3.0.2 in PostgreSQL 9.6.


texts: Using planet.bz2 with ArcGIS Editor for OSM, sync'ing changes

texts:
Overview: We are attempting to clean up and revise the OSM data of our County by comparing it to a more accurate (public domain licensed) street layer within the ArcGIS Editor for OpenStreetMap.  First, we downloaded a planet dump .bz2 file of California roads from Cloudmade.  Next, using Osmosis and the approximate coordinates of our County, we clipped the .bz2 to a more manageable size for quicker loading into ArcGIS.  Finally, using OSM Simple Loader tool, we brought the .bz2 into a feature dataset in ArcGIS.

Question: How can we properly set up Geodatabase tables for the feature dataset to be able to synchronize the .bz2 dump with the most current data on OSM servers? I understand that it requires a specific schema to work properly, we want to make sure we can update the planet dump, so we are editing against the most recent data.  We need help with this final step.


texts:
How can I configure base layers in Geonetwork 3.0.2?

I have tried to edit config-gui.xml and config-viewer.xml, but nothing happens.

Does it only work for the old version?


texts: ArcPy and File Geodatabase Character Encoding

texts:
I have a JSON file containing special characters, which can be represented in Unicode. However, as soon as I export that JSON file to a file geodatabase using arcpy, the special characters all "break," i.e., they come out as something different. A Google search tells me (via this table listing string representations of special characters) that the character I'm having a problem with in my sample dataset is in UTF-8/Latin-1.

The code I'm using for this is:

arcpy.JSONToFeatures_conversion(jsonFileName, os.path.join(fgdbPath, outputName))

I have also specified the encoding for the Python (2.7) file via:

import sys
reload(sys)
sys.setdefaultencoding('utf-8')

I have also varied the above with "latin-1" and various other specifications for the encoding, mostly experimentally as I have little experience with encoding.

But regardless, where the content of my JSON file is for example "Búzi", the content of my feature (exported to an FGDB) will be "BÃºzi".

Nowhere can I find out about a parameter to pass to JSONToFeatures_conversion, or to CreateFileGDB_management (where I create the FGDB in question), or an Arc environment variable to set, which will allow me to specify the desired encoding.

How do I export to FGDB with the desired encoding (or providing the desired result by whatever means)?


texts: Subtracting image value from corresponding month from image collection in Google Earth Engine

texts:
I have extracted sea surface temperature(SST) images from Google Earth Engine for a specific time period and region
var start = ee.Date(2012-01-01);
  var end = ee.Date(start).advance(1,'month');
  var date_range = ee.DateRange(start,end);

var sstImg = ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI')
.filterBounds(reefs_500) //Filter by Study Area
.filterDate(date_range) // startdate, enddate
.select('sst')
.sort('system:time_start', false);

I want to calculate the long term mean SST per month for the region and then subtract that value from the image collection above.
I have calculated the mean for each month from 2002-2012 and have this as an image collection for 12 different images one per month.
var sst2002_12 = ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI')
    .select('sst')
    .filterDate(ee.Date('2002-01-01'), ee.Date('2012-12-31'))

var months = ee.List.sequence(1, 12);

var byMonth = ee.ImageCollection.fromImages(
      months.map(function (m) {
        return sst2002_12.filter(ee.Filter.calendarRange(m, m, 'month'))
                    .select(0).mean()
                    .set('month', m);
}));

How can I subtract only the mean from the image for the same month e.g. should the image come from January, it takes the Jan (or 01) value from the image collection?
I have this code to subtract a specific value from each image.
var series_daily = sstImg.map(function(image) {
    return image.subtract(mean)
    .set('system:time_start', image.get('system:time_start'));
});

But not sure how to adapt it so it only picks the correct value to do the subtracting from the list.


texts: Where to download CSW Client for ArcMap (ArcGIS for Desktop)?

texts:
I just tried to find and download the CSW Client for ArcMap so that I could use it with ArcGIS 10.3 for Desktop to test an answer at Is there public OGC Catalog Service (CSW) available which lets ArcGIS for Desktop with CSW Client Add (WMS) To Map?

I located a Geoportal Server Downloads page on GitHub and clicked on the Desktop Clients for ArcGIS 10.2 link because there was not a link to Desktop Clients for ArcGIS 10.3.

When I tried to install it all the same I was blocked by a message saying:

  CSW Client for ArcGIS Desktop requires ArcGIS 10.2 for Desktop CSW

  Client for ArcGIS Desktop will not be installed at this time.

Is there a CSW Client for ArcGIS 10.3 Desktop available for download anywhere?


texts:
I am new to image processing and coding and I am having trouble applying a lowpass and highpass filter to a NAIP image. I have practiced with Landsat 8 images with a successful output but I am getting a white/ blank image for the NAIP image. I know I might be overlooking some minute detail, but any help is appreciated. I tried searching for more GEE documentation on this matter but found very limited information.
Here is a sample of my code:

var image=ee.Image('USDA/NAIP/DOQQ/m_3110613_se_13_060_20190226_20190504').select('R','G','B');
// Region of interest for section with buffer of 2.5 km
var roi = ee.Geometry.Point(-106.4126, 31.7835).buffer(2500);

Map.setCenter(-106.4126, 31.7835,15);
Map.addLayer(image,{},'RGB NAIP Image: Q3');

// lowpass filter parameters
// Define the low-pass kernel.
// This kernel creates a 7X7 square kernel.
var lowpass = ee.Kernel.square({radius: 7, units: 'pixels', normalize: true});

//Smooth the image by convolving with the low-pass kernel.
var smooth = image.convolve(lowpass);

//Display the filtered image
Map.addLayer(smooth, {bands:['R','G','B'], max: 0.5}, 'smoothed');


texts: Temporal controller - time-dependent symbology based on @map_start_time

texts:
I was just recently whatching Nyall Dawson's video "QGIS 3.14 time handling and map animation demo".
at min 47 (nyall dawson / youtube video min 47 ) things are getting interesting from my point of view.
While testing the "temporal controller" features I've had some confusing results when trying to implement the time dependent symbology approach shown in the video - imag_01 shows my results following the video-tutorial:

data: a GPS track from a yacht race in the Adriatics. data
acquisition interval is about 5 min
point symbols (yellow): in the foreground are controlled by QGIS'
temporal controller
points (grey): in the background are non-temporal and displayed just
for giving an overview
temporal controller settings are:
step width: 15min (= 0,25 hours)
the displayed interval starts at 2018-10-15 15:45:00 and ends at 2018-10-15 16:00:00 (see the values next to "Frame:" entry)

point size is depending on a temporal expression - just like in the video:
hour(@map_start_time - time_s_oesec") (see the expression in image_02)
where

@map_start_time is the QGIS-internal variable - which should be equal to the Frame start value, i.e. '2018-10-15 15:45:00' in my example
time_s_oesec is a reference to the temporal field in my data - i.e.
the location of the gps-tracking timestamp data
(... values in the field time_s_oesec are stored without seconds in this example to keep things simple)

Since @map_start_time is the lower limit of the displayed time interval any timestamp within the interval is greater or equal @map_start_time. Scaling the yellow points (according to Nyall Dawson's video) is based on the above-mentioned expression: hour(@map_start_time - time_s_oesec).
image_01 shows my test-labels

1st line:  GPS point field value e.g. "time_s_osec:        15.10.18 15:50:00"
2nd line: current value of the map start variable e.g. "@map_start_time:   15.10.18 15:45:00"
3rd line:  "@map_start_t* - time_s*:     1,92"
that values should be negative due to the expresion @map_start_time - time_s_oesec
and they should NOT range from > 1 to > 2 within the defined interval
4.line: this is my check-value --> this value is the output from:
hour(to_datetime('2018-10-15T15:00:00') - "time_s_osec" )
where '2018-10-15T15:45:00' is a fixed text which is transformed to hours based on QGIS datetime-format and should exactly correspond to the actual value of the variable @map_start_time

Obviously the expression in line 4 returns valid results - see/compare the following examples in image_01:
For the point at gps-time in the map 2018-10-15T15:45:00

to_datetime('2018-10-15T15:45:00') /  time_s_oesec: 15:45  / result =
-0,00  ... which is as expected  but @map_start_time - time_s_oesec returns a value of 2,0

for the point at gps-time in the map 2018-10-15T15:29:00

to_datetime('2018-10-15T15:45:00') / time_s_oesec: 15:29 / result =
0,27  but @map_start_time - time_s_oesec returns a value of 2,27

Summary:  expressions in line 4 - which are based on the text-check-value '2018-10-15T15:45:00' return positive values for gps-timestamps before "frame start time" 15:45 (since the GPS timestamp is less than the start time) and negative values after 15:45 since it is the other way round.
While the expression hour(to_datetime('2018-10-15T15:00:00') - "time_s_osec" thus results to correct values the way the "temporal manager" handles @map_start_time in that context is strange/faulty.
But... Is it worth a bug report ?
If needed I can provide test-data
image_01:

image_02:


texts:
I have a problem when using the con() argument.

What I want to do is to remap all raster of a raster file ("Red_SL") cells with values larger than 40 to the new value 5. All other values should be replaced by another raster file called "remap2". "remap" 2 contains only 4 definite values (1,2,3,4), and "Red SL" has continuous values between 0 and 90.

My syntax is the following:

import arcpy
from arcpy.sa import *

arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"D:\\Projekt\\Testprojekt\\Daten\\Ergebnisse\\"
arcpy.CheckOutExtension("Spatial")
remaped2 = "remaped2.img"
slope = "Red_SL"

CF= Con((slope < 40), 5, remaped2)
CF.save("Red_SL_Neigung_int.tif")
arcpy.CheckInExtension("Spatial")

now the problem is: I do not get an error message but the produced file is nonsense!
Python is somehow not able to recognize the slope value in: CF= Con((slope < 40), 5, total)
This means it doesent matter what value I give python (here I set it to 40, but I will get the exact same result if I give values of 1 or 3000), the result will always be the a new file which is exactly the same as "remaped2". Python therefore thinks all values are below 40 (or whatever number I give), which is not true.

The exact opposite happens if I change the "<" to ">" (CF= Con((slope > 40), 5, total)).
Now I will get i new file, which has only 5 in every raster cell. Here again the value whch I type in (instead of 40) does not matter, any value between 1 and 3000 will produce the exact same result. I double checked, but the file "Red_SL" has all values between 0 and 90, therfore the result I get from python is nonsense.

I already tried (without sucess):
-Change "Red_SL" from floating to integer
-try differend formats (tif, grid, img) for "Red_SL"
- replace comma in floating raster with points (german vs. englisch delimiters)

Can anyone see what is wrong with my use of the Con arguement?


texts: Using REST API to get the data collected in Survey 123?

texts:
I'm trying to figure out how I can get the data collected with Survey 123 out of ArcGIS Online via the REST API.

I am using Postman to try and formulate the POST request, but, all it shows me is the feature service metadata. In fact, it returns JSON that contains everything about the feature service but no actual data.

What am I doing wrong?


texts: Remove duplicate names in FeaturePopups cluster tooltips and popups

texts:
I am using FeaturePopups control to show tooltips and popups for clusters.
I have features connected (ends of road segments) with the same names.
Is it possible to remove them from tooltips and popups, so to have a list with unique names?

Thank you,
Bogdan


texts: QGIS: \ in my expression gets returned as a?

texts:
Right now I am building a expression in a column I am gonna be using later. You see I have a problem where a expression in one of my columns this one:
'\gpkg_content.dbf'

Gets returned as this.
?pkg_content.dbf

So my \ is always converted to a ? while clearly I did not clarify that.
Why does this happen?


texts: Creating standalone application with same directories and qgis-auth.db

texts:
I try to have a standalone QgsApplication (api documentation), but, instead of a lot of code samples found like this one below (I'm on Ubuntu):

from qgis.core import QgsApplication

QgsApplication.setPrefixPath("/usr", True)
app = QgsApplication([], False)
app.initQgis()

I want the same result for print(app.showSettings()) as in the QGIS Desktop application.

In QGIS Desktop :

>>> print(QgsApplication.instance().showSettings())

État de l'application
variable d'environnement QGIS_PREFIX_PATH :
Préfixe:        /usr
Chemin des extensions:      /usr/lib/qgis/plugins
Chemin des paquets de données:  /usr/share/qgis
Nom du thème actif: default
Chemin du thème actif:  /usr/share/qgis/resources/themes/default/icons/
Chemin du thème par défaut: :/images/themes/default/
Chemins de recherche SVG:   /usr/share/qgis/svg/
        /home/me/.local/share/QGIS/QGIS3/profiles/default/svg/
Chemin de la BD utilisateur:    /usr/share/qgis/resources/qgis.db
Chemin de la BD d'authentification: /home/me/.local/share/QGIS/QGIS3/profiles/default/qgis-auth.db

In standalone Python script :

>>> from qgis.core import QgsApplication

>>> QgsApplication.setPrefixPath("/usr", True)
>>> QgsApplication.setThemeName("default")
>>> QgsApplication.setAuthDatabaseDirPath("/home/me/.local/share/QGIS/QGIS3/profiles/default/qgis-auth.db")
>>> app = QgsApplication([], False)
>>> app.initQgis()
>>> print(app.showSettings())

Application state:
QGIS_PREFIX_PATH env var:
Prefix:     /usr
Plugin Path:        /usr/lib/qgis/plugins
Package Data Path:  /usr/share/qgis
Active Theme Name:  default
Active Theme Path:  /usr/share/qgis/resources/themes/default/icons/
Default Theme Path: :/images/themes/default/
SVG Search Paths:   /usr/share/qgis/svg/
        /home/me/.local/share/profiles/default/svg/
User DB Path:   /usr/share/qgis/resources/qgis.db
Auth DB Path:   /home/me/.local/share/profiles/default/qgis-auth.db

It creates a new profiles directory and it's not what I want.

I also tried to initialize the db auth Manager with no results, it leads to a Python segmentation fault.

How to have the same Auth DB Path in QGIS Desktop and in a standalone Python script?

In summary, how to have (done, see the edit) :

Auth DB Path:  /home/me/.local/share/QGIS/QGIS3/profiles/default/qgis-auth.db

Working on Ubuntu 18.04, QGIS 3.8.

With the proposition of @GermánCarrillo:

call QgsApplication.setAuthDatabaseDirPath("/home/me/.local/share/QGIS/QGIS3/profiles/default/") (without qgis-auth.db) after app.initQgis()

it works but it failed after :

>>> from qgis.core import QgsApplication

>>> QgsApplication.setPrefixPath("/usr", True)
>>> QgsApplication.setThemeName("default")
>>> app = QgsApplication([], False)
>>> app.initQgis()
>>> QgsApplication.setAuthDatabaseDirPath("/home/me/.local/share/QGIS/QGIS3/profiles/default/")
>>> print(app.showSettings())

[... same as before ...]
Auth DB Path:   /home/me/.local/share/QGIS/QGIS3/profiles/default/qgis-auth.db

Ok, I have the correct Auth DB Path but if I try to get the QGIS app authManager, the database path is still not correct :

>>> authMgr = app.authManager()
>>> authMgr.authenticationDatabasePath()

'/home/me/.local/share/profiles/default/qgis-auth.db'

I thought in the first place that if the Auth DB Path of the QGIS app is properly configured, the database path of the app authManager would be too.

Now, how to have :

app.authManager().authenticationDatabasePath() equals to /home/me/.local/share/QGIS/QGIS3/profiles/default/qgis-auth.db

or how to get access to the Desktop app credentials / authentication database ?


texts: Copying all coordinates at once in QGIS?

texts:
I had a line with points and I want to get all of its coordinates by once,
The only way that I can copy the coordinates by this step
(Identify feature> right click the point that I choiced ,copy feature and paste it to get the coordinates,

but this step only selecting a coordinate , i want to get all coordinate,


texts: I have a city road network. I have to divide sections that are separated by intersections into one level. With ArcGIS this is done with Planarize. I tried split commands, but it applies to all sections of the net. How can I do this by plotting specific sections?


texts: How do I display NYC community specific data on an interactive map without longitude/latitude coordinates?

texts:
I have downloaded the NY Department of Sanitation's Monthly Tonnage Collection data from https://data.cityofnewyork.us/City-Government/DSNY-Monthly-Tonnage-Data/ebb7-mvp5#revert and I want to use this information to create an interactive map using CartoDB.

The problem is this data set does not contain long/lat coordinate data and thus cannot be mapped. Also, the smallest spatial denomination provided in the data set is community/district which means I probably need several coordinates to map each community/district.

I'm new to making interactive maps and GIS in general. What's the best way for me to proceed? Should I pluck out appropriate coordinates for each community/district and push them into the data set?


texts: A line between point to nearest neighbour point using QGIS

texts:
I have a layer of points in QGIS which represent electric poles

How do I create a polyline layer that would have lines running between the poles so it can represent the the electric cables that connecting the poles together?

I have been searching on the web and could not find an answer.


texts: Recommended way to explicitly distinguish 2D from 3D locations

texts:
In strongly typed object-oriented programming environment, I want to make an explicit distinction between two types of Locations:

One type would be 2D (along Earth surface), thus having Latitude and Longitude. With a DEM available, Elevation would be a function of Latitude, Longitude, and the given DEM model. This type would represent positions of terrestrial or water-floating objects moving along the Earth surface.
Other type would be (inherently) 3D, with independend Elevation, intended to represent aerial, underwater or subterraneal positions and trajectories.

So my question is:

  Is there already any de-facto standard or commonly used name for such distinct types of locations?


texts: Running Fishers (1996) 'Local Offset Viewshed' with ArcGIS Desktop?

texts:
Can ArcGIS Desktop run the Local offset viewshed as described by (Fisher, 1996) on page 1300-1301 below?

There's a similar operation called 'Invisibility depth' available through the plug-in 'Advanced viewshed analysis' for QGIS, but I'm looking for an alternative for comparison purposes.

Fisher, Peter F. "Extending the applicability of viewsheds in landscape planning." Photogrammetric Engineering and Remote Sensing 62.11 (1996): 1297-1302. http://info.asprs.org/publications/pers/96journal/november/1996_nov_1297-1302.pdf


texts: Iterating extraction of polylines and annotations from CAD to Esri Geodatabase

texts:
I have successfully converted a DWG lines and annotation in a geodatabase annotation feature and polyline feature using Arcgis 10.2.2
And now I want to construct a model builder to automate and iterate this process for my 19000 remaining dwg.
I started to build it but I am confuse with which tools to use so that the model works and names of feature doesn't get duplicated inside geodatabase.

The steps I used to convert it manually one by one is:

1-Create a file geodatabase

2-Open ArcMap and insert the annotation and the polyline layer of the Autocad file in Arcmap

3-Click on annotation layer and convert to geodatabase annotation

4-Click on polyline layer and export data to the geodatabase

5-Then from toolbox merge all the converted annotation layer in one annotation layer and do the same for the polyline layer

Now I have created a geodatabase that I want to use for this conversation and I started building my model builder but I am very confuse with which tool I can automate it, suggestion change something even to the one I already started because I am not sure if what I started is correct? Below is the screenshot of my model builder


texts: Is GeoExt styler extension available for Geoserver 2.3.3?

texts: I've switched from my older version of geoserver to geoserver 2.3.3. In the older version, the  styler worked just fine, but I could'nt find a styler extension that matched the new version. I tried with styler of 2.1.7 version in geoserver, but it failed to work. Can anyone please tell me how to get a styler extension that works with the new version. Kindly provide download links if available.


texts: Run script for many profile path at once in global mapper

texts:
This code when I run it in Global Mapper, that create a .xyz file of profile path between two points

GENERATE_PATH_PROFILE \
FILENAME="C:\Users\Bella\Desktop\ProfPath.xyz" \
ELEV_UNITS=METERS \
POINT_SPACING=10.0 \
START_POS=765287.9,3813969.2 \
END_POS=771481.8,3815738.6 \
SAVE_DIST_Z_FILE=YES \

1- If I would to create many segments at once, all the segments share the same origin point, how to do it in script, so I can export it all in a file?

2- Could I use the geographic coordinates in START_POS and END_POS instead of utm projection?


texts: gdal_translate tiles won't save anything

texts:
Other Tiles question

I copied the code and tried to run it, it runs the for loop the cmd says "100.. done" but there aren't any files in the set directory. Could be there any problem with the permission? If yes, how can I elevate a QGIS Python addon?

 def start_splitting_tiles(self):
    in_path = self.dlg.currently_choosen.toPlainText()
    input_filename = in_path

    while (input_filename.find("/") != -1):
        input_filename = input_filename[input_filename.find('/') + 1:]

    print "input: "+input_filename

    out_path = self.dlg.outpath.toPlainText()

    #if(out_path.find('\\') != -1):
    #   out_path = out_path.replace("\\", "/")

    out_filename = 'tile_'
    print "output: "+out_path

    tile_size_x = 1500
    tile_size_y = 1500

    ds = gdal.Open(in_path)
    band = ds.GetRasterBand(1)
    xsize = band.XSize
    ysize = band.YSize

    for i in range(0, xsize, tile_size_x):
        for j in range(0, ysize, tile_size_y):
            com_string = "gdal_translate -of GTIFF -srcwin " + str(i) + ", " + str(j) + ", " + str(
                tile_size_x) + ", " + str(tile_size_y) + " " + str(in_path) + " " + str(
                out_path) + str(out_filename) + str(i) + "_" + str(j) + ".tif"
            os.system(com_string)
            #os.system(input())


texts: GRASS 6.4.4 vector command in crashes QGIS 2.4 (g.gisenv.exe has stopped working)

texts:
I am running Windows 8.1 64 bit, and have QGIS 2.4 and 2.6 installed. When I try to use a GRASS 6.4.4 vector command (specifically v.clean), I receive an error message in an endless loop, and the only way to cancel it is to close QGIS. This happens for both 2.4 and 2.6, and I have tried uninstalling and re-installing GRASS.

This is the error message I recieve:

and this is how I have GRASS set up:


texts: Creating pre-determined size polygon with average raster value using ArcGIS Desktop?

texts:
I have a simple raster with multiple pixels (5 million+), and each pixel has a value from 0-1 representing % values. So one pixel might be 0.361 (36%) and another might be 0.053 (5%).

Is it possible for GIS with something like a moving window to create polygons for me that have a specified value from the raster?

I wanted 3 polygons of a certain size (acres) where each has an average value of 70%. Each pixel in the raster is 1 meter.

I could do this manually in an iterative way by drawing polygons than calculating statistics, then re-drawing, but can GIS do it for me?

I am not technically familiar enough with GIS to know, and google has indicated that I might need to go about this in a round-about way.


texts:
I must apologize in advance, I'm just getting into GIS and trying to figure things out.

I have an SQLite spatial database that I'm pretty sure was created in TatukGIS. I've tried to open it in QGIS but the database format is different than that of a Spatialite database. I would like to find a way to convert the SQLite database created in TatukGIS to a Spatialite Database that can be opened in QGIS.

I know it has something to do with the OpenGIS simple features specification and how Spatialite may not fully comply with those specifications, but I don't have enough experience to know how to go about converting it (if it's even possible).

Here is a link to the Databases to give you an idea how they are organized.

https://drive.google.com/folderview?id=0B2TMeZ6iEFvHa3BBOVJpVVpBY1E&usp=sharing

TatukGIS.sqlite is in epsg:3569,
Spatialite is in epsg:3857


texts: Landsat 8 to the temperature map (in ArcGIS 10.0 -10.2.2 or ENVI 5.0)

texts:
I need to create a temperature map from Landsat 8. I found  an algorithm that explains how to do it in ArcGIS:

http://blogs.esri.com/esri/arcgis/2014/01/06/deriving-temperature-from-landsat-8-thermal-bands-tirs/

Unfortunately I can’t use this method, because the “local function” (that is used in Step 3  - “To calculate the natural log, you need to use local function”) is blocked in my version of ArcGIS. The program gives: NOT LICENSED "This raster function requires a license"
It’s strange, because now I use the trial version of ArcGIS 10.2.2 (that gives all possibilities for 2 months), maybe I can find this function somewhere?

Also I tried to realize this algorithm on the permanent version of ArcGIS 10.0 But here I can’t find the Raster Function Template Editor" (from Step 1 – B) (I don’t have any “Commands” in Customize >>Raster).
So, two problems in the both versions of ArcGIS  that I can’t solve.

But as I understood, this method explains how to create a Function chain which can be used for any Landsat 8 for preparing a temperature map.  Maybe some special script for ArcGIS for  this purpose  has been developed?

Or maybe you can advise how I can create a temperature map from Landsat 8 in ArcGIS (by dint of this method or other).
Also I have ENVI 5.0,  so I’ll  be glad advises for method by dint of ENVI as well.


texts: Is there a way to autocorrect miss-spelled cities for geocoding?

texts:
I could not geocode a subset of my dataset (around 140,000 observations) which contains full postal addresses from Germany, i.e. zip code, city and street. Apparently, the geocoder failed to convert the addresses due to mistakes in the data. Eyeballing led me to the impression that the most prevalent reasons are the following

minor spelling mistakes in terms of the city name. Example: "Neunburg v. Wald" instead of "Neunburg vorm Wald" or "Hessheim" instead of "Heßheim"
missing part of the name. Example: "Pohlheim" instead of "Pohlheim-Watzenborn"
wrong match between zip code and city name (maybe because the address was recorded before a change in the zip code occured, e.g. two zip codes were merged)

This is why I would like to "auto-complete" / "auto-correct" the city names. I could imagine to do the following: If I had access to a database that contains all zip codes and cities for Germany, I would make a list of "suggestions" among which I choose the one that is "closest" to the wrong city name. The way one could go about the list of suggestions is to match all cities that have the same plz (not a one-to-one mapping) or match those cities from the database that start with the same first x characters (under the assumption that the error does not occur in that range). Then, one could pick the closest string (city name) based on an approximate string match algorithm like the Levenshtein algorithm.

This leads me to two subquestions:

Is there a possiblity to extract this information (zip codes and city names) from the OSM database (I do not have an installed instance though)?
Should I rather rely on service providers like Google?


texts: How to perform a spatial join of point and polygon layers in PostGIS?

texts:
I have a point dataset representing households that I want to associate with a parcel layer (i.e. assign its parcel_id). With ArcMap, I could spatially join the polygons to the points and specify that the join have a certain search radius and use the nearest polygon.

I wish to perform this with postgres and was wondering what the most prudent way of doing this was. I can't use ST_contains because points fall on streets and not on actual parcels. In stead of ST_contains, is there another function allowing for proximity checks?

update households h set parcel_id = p.parcel_id
from parcels p where ST_contains(p.the_geom, h.xy);


texts: Translating EPSG:4326 lon lat coordinates to xy using pyproj

texts:
I'm trying to work out the details of using pyproj to translate the lat lon coordinates from the USGS Earthquake feed to xy coordinates corresponding to those positions on a basemap of a determined width and height and origin at the upper left. My code so far naively uses PIL/pillow to translate lon lat to xy positions and draw circles without taking into account the EPSG:4326 projection of the basemap:
from PIL import Image
from PIL import ImageDraw

# open basemap image file
basemap = Image.open(basemap_path).convert('RGBA)

# resize to desired map size
basemap.thumbnail(width, height, Image.LANCZOS)

# get proportional height
width_bmp, height_bmp = basemap.size

# create background frame and paste basemap on it
img=Image.new('RGB',(width, height), color = '#000000')
img.paste(basemap, (0,0), basemap)
draw = ImageDraw.Draw(img, 'RGBA')

width_scale = width/360
height_scale = height_bmp/180

# usgs data has been parsed into a list
for quake in earthquake_list:
    lon = float(quake["longitude"])
    lat = float(quake["latitude"])
    mag = float(quake["mag"])

    # want to use pyproj to translate coordinates here instead of the following
    cx = (lon + 180) * width_scale
    cy = (90-lat) * height_scale
    r  = scaleRadius(mag)

    # draw earthquake circles
    draw.ellipse((cx-r, cy-r, cx+r, cy+r), fill = colormap(mag))

draw = ImageDraw.Draw(img)
img.save(filepath, quality=100)

I've gone through the documentation and have a general sense that I would use something like:
cx, cy = pyproj.transform("EPSG:4326", "xy", lon, lat)

with some width and height dimensions. Running the line above gives the error:
TypeError: p1 must be a Proj class

Update I do get values for:
p = Proj(proj = 'longlat', ellps='WGS84')
cx,cy = p(lon, lat)

In this case a (lon, lat) input of (-110, 39) gives the output 1.923, 0.688 which isn't right. What I need is a projection function that takes the lon, lat and the width and height of the image and outputs the x, y coordinates. I would write this on my own, but I can't find a suitable equation for that conversion and was hoping to find this capability in proj.
UPDATE Here is the basemap I am using.

And here is the resulting drawing. I'm assuming the centerpoints for the circles marking earthquakes are off, because I am using Cartesian coordinates for lon, lat, but the basemap was created with a projection and thus must have a distortion from x, y Cartesian coordinates.

Also, I realize I need to put my legend on a logarithmic scale. Haven't gotten to that yet.


texts: Read a table from an ESRI file geodatabase (.gdb) using R

texts:
I am trying to read a table directly from an ESRI file geodatabase into R.  An example data file can be downloaded here.  The database contains a point feature class (Zone9_2014_01_Broadcast) and two linked tables (Zone9_2014_01_Vessel and Zone9_2014_01_Voyage).  You can read the shapefile in R using readOGR from the rgeos package:

library(rgeos)
library(downloader)

download("https://coast.noaa.gov/htdata/CMSP/AISDataHandler/2014/01/Zone9_2014_01.zip", dest="Zone9_2014_01.zip", mode="wb")
unzip("Zone9_2014_01.zip", exdir = ".")

#  Not Run (loads large point file)
#  broadcast <- readOGR(dsn = "Zone9_2014_01.gdb", layer = "Zone9_2014_01_Broadcast")

The two linked tables also show when you use ogrListLayers or ogrInfo.  However, ogrInfo gives a warning:

  Warning message:
  In ogrInfo("Zone9_2014_01.gdb", layer =  "Zone9_2014_01_Vessel") :
    ogrInfo: all features NULL

And if you try to use readOGR on the tables you get an error:

vessel <- readOGR(dsn = "Zone9_2014_01.gdb", layer = "Zone9_2014_01_Vessel")

  Error in readOGR(dsn = "Zone9_2014_01.gdb", layer =  "Zone9_2014_01_Vessel") :
    no features found
  In addition: Warning message:
  In ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv,  :
    ogrInfo: all features NULL

Thus, it appears that only geographical features can be read by readOGR.  Is there any way to import the tables directly into R or is the only solution to first export them from ArcGIS as *.dbf (or *.txt) files as in this answer?

As an addition, if anyone can provide calls from R to a python script that automates the export of *csv (preferrably) or *.dbf files, that would be an acceptable work-around.  The solution just needs to be scalable and automated.


texts: arcpy.NumPyArrayToFeatureClass exception: "TypeError: narray.fields require"

texts:
I have a database application that I connect to using SQLAlchemy. In the database, I store coordinates in the fields x_coord and y_coord, with a third field spatial_reference_code that stores the EPSG code of the coordinate system that x_coord and y_coord are in for the record.

I've written a generic function that takes a SQLAlchemy query object with any preexisting filters on the data table and exports it to a feature class. The function converts the SQLAlchemy object to a pandas data frame, then uses the pandas data frame (after using data_frame.as_matrix()) as the data  parameter to arcpy.da.NumPyArrayToFeatureClass. When running this code, it fails in NumPyArrayToFeatureClass with the error "TypeError: narray.fields require"

Below is a minimum reproducible example, adapted from Esri's sample code, that generates the problem:

import os
import arcpy
import numpy
import pandas

out_fc = os.path.join(os.getcwd(), "test.shp")

# Create a numpy array with an id field, and coordinate fields

array = numpy.array([(1, -116.809, 34.271),
                     (2, -114.241, 33.684)],
                    numpy.dtype([('idfield',numpy.int32),('x_coord', '<f8'), ('y_coord', '<f8')]))

# converting to pandas data frame to simulate conditions in my actual code - this causes the failure
df = pandas.DataFrame(array)

# Define a spatial reference for the output feature class
sr = arcpy.SpatialReference(4269)

# Export the numpy array to a feature class using the X and Y fields to
#  represent the output point feature
arcpy.da.NumPyArrayToFeatureClass(df.as_matrix(), out_fc, ['x_coord', 'y_coord'], sr)

My code does not start with a numpy array and convert it to a data frame before converting it back, but the example does that to simulate the conditions of my code (where the data is retrieved from SQLAlchemy into a pandas data frame). If I run that sample code and pass "array" directly to NumPyArrayToFeatureClass, it works - but when I convert it to a pandas data frame and then back using df.as_matrix() - it fails.

For machine specs, this code is running in the Python interpreter provided by ArcGIS Pro 1.3 (Python 3.4.4 - in a conda environment I created from the main ArcGIS Pro environment). Pandas is version 0.17.1, numpy is version 1.9.3, and SQLAlchemy is version 1.1.2. Running on Windows 10 (10.0.14393).

I see one other reference to this issue over on StackOverflow, but that post does not resolve my problem.

Has anyone worked with this arcpy.NumPyArrayToFeatureClass and can help me debug this issue?


texts: MapInfo cannot join after deleted some rows with "del"

texts:
I tried to join my table with another one, with the update column. But it gives me an error, that it cannot be updated. Actually it does update the table to a certain point.

I deleted some rows before I tried to update the table, manually by clickling left on the row and pressed the key "del", so the object got removed, the row still kinda exist, but only in grey. But that is the point where the join stops, from there on there are no further entries.

I tried to close and reopen the table, but the grey rows are still in there, any idea?


texts: All (GeoTIFF) raster values NA after converting to velox

texts:
Cross-posted to the velox GitHub repo:  https://github.com/hunzikp/velox/issues/40

I'm able to extract average tree cover values for each of the three small counties in the US State of Delaware, but I don't seem to have enough CPU or RAM to do the same thing for all counties in the US.

I'm getting my tree cover data from https://prd-tnm.s3.amazonaws.com/StagedProducts/Small-scale/data/LandCover/tree48i0100a.tif_nt00840.tar.gz.  See also https://www.sciencebase.gov/catalog/item/581d0548e4b08da350d52653

library(raster)
library(tigris)

tigris.state.counties <- counties(state = 10)
tree.raster <-
  raster("tree48i0100a.tif")

polygon.trees <-
  extract(tree.raster,
          tigris.state.counties,
          fun = mean,
          na.rm = TRUE)

report <-
  cbind.data.frame(tigris.state.counties$NAME, polygon.trees)

print(report)

  tigris.state.counties$NAME polygon.trees
1                 New Castle      70.75920
2                     Sussex      68.17089
3                       Kent      67.00225

So I have set up an R environment on a 64 GB AWS instance, and I'm trying to take advantage of the velox fast extraction.  But all of the values I extract are NA.

What am I doing wrong?

UPDATE:  I still get NAs even if I use the small polygon option

FWIW, I was able to replicate the velox extract tutorial and even modify it to get the night light intensity in the Delaware counties.

My attempt to extract tree cover in Delaware counties with velox:

library(velox)
library(tigris)

tree.velox <- velox("tree48i0100a.tif")

tigris.state.counties <- counties(state = 10)

tree.mean.mat <-
  tree.velox$extract(
    sp = tigris.state.counties,
    fun = function(x)
      mean(x, na.rm = TRUE)
  )

report <-
  cbind.data.frame(tigris.state.counties$NAME, tree.mean.mat)

print(report)

     tigris.state.counties$NAME tree.mean.mat
591                  New Castle            NA
1150                     Sussex            NA
2861                       Kent            NA


texts:
I'm trying to symbolize a line layer in CartoDB using CartoCSS. I need my line to be a series of diagonal lines, something like this -> /////////////

In QGIS I would do this as a marker line with a straight line set at an angle for the marker. However, I can't figure out whether such functionality is present in CartoCSS. I can't see anything in the documentation for marker lines. Is this not possible?


texts: The shapefiles of two different contours overlap with each other

texts:
I tried to create two contour lines with a relatively small difference (0.5meter) on a DEM dataset.
For example, I created an 8-meter contour(the red one) and then created an 8.5-meter(the blue one) contour. Their shapefiles overlapped with each other at the place which should be the position of 8-meter contour.
What's weirder is that I cannot select the 8.5-meter shapefile while I can select an invisible line of 8.5-meter beside its shapefile (highlighted blue line).

Here's the code I used to create these contours, I set the interval as 200 so it will only draw one line:

import arcpy

from arcpy import env

env.workspace = "D:\Arc GIS\"

arcpy.Contour_3d("Dem.tif","D:\Arc GIS\Dem,-8.shp", 200, -8)

arcpy.Contour_3d("Dem.tif","D:\Arc GIS\Dem,-8.5.shp",200, -8.5)


texts: How can I perform a Join by Location for a large vector and very large point file?

texts:
I have no programming skills, but am thinking I may need to at least figure out how to run this one script from inside QGIS.

I have a vector file of global quarter degree grid cells and a point file containing an incredibly large number of points with population values.

I want to create a polygon file that has the sum of the population values (not just the count of the points) for each grid cell.

This is really just a basic Join by Location (Spatial Join), but the files are so large it will not work.

Can anyone help me?


texts:
Can the following be done in ArcMap? It should be formula based without manual formatting as performed by my predecessor (in MapInfo using a lot of manual steps) to get the following output.

So

Grave_No on right parallel to the y axis of the grave.
Name on top center parallel to the x axis of the grave.
Date of Death on bottom center parallel to the x axis of the grave.

I expect that we will need to copy the layer three times to apply custom labels for each of these situations

In MI
-I can't autorotate the labels and place them on the left/make sure the angle is the same as the grave. I know it can be done manually per label but this is far to much work.

In QGIS
- See Custom labeling with same angle as rectangular object and specific position using QGIS?

In ArcMap
- In maplex, it seems as if we can only rotate point labels but I have a polygon shape. We can't use centroids as it will go out of the grave site.

So QGIS seems to get the closest...and while I am happy to have a QGIS solution I have asked a new question (as there were no responses) to just have an ArcMap solution.

Note 1. The name sometimes does not show but the Date of Death does, in other cases this is reversed. Both should be shown. 2. Some will have 2 or more people in the grave so the data will be there twice - so grave no will be E-062a and E-062b for example with 2 names and dates of death as well.

So can I do all the above in ArcMap? the example shows how far I could get it.


texts:
I have used OpenLayers2 in order to show a layer with points from Geoserver.
I try now to transform the code from Openlayers2 to OpenLayers3 but I am facing some issues.
I have walked through this nice tutorial: OpenLayers3
and I have a good understanding of how to add layers from localhost geoserver in OL2.

This is part of my code:

function init() {

map = new ol.Map({
    target:'map',
    renderer:'canvas',
    view: new ol.View({
        projection: 'EPSG:900913',
        center:[0,0],
        zoom:5
    })
});

var newLayer = new ol.layer.Tile({
    source: new ol.source.OSM()
});

map.addLayer(newLayer);

var vectorLayer = new ol.layer.Tile({
source: new ol.source.TileWMS({
        preload: Infinity,
        url: 'http://localhost:8080/geoserver/gwc/service/wms',
        serverType:'geoserver',
        params:{
            'LAYERS':"dSpatialAnalysis:categoriesdata", 'TILED':true
        }
    })
});
map.addLayer(vectorLayer);
}

For some reason I get a:

GET http://localhost:8080/geoserver/gwc/service/wms?SERVICE=WMS&VERSION=1.3.0&R…TYLES=&BBOX=2504688.542848654%2C0%2C3757032.814272982%2C1252344.2714243277 400 (Bad Request)

I understand that somethings regarding the projections have changed in OL3 and I wonder if this is related with the errors I get.
I also checked if the layer works correctly. I can view it in the "Layer Preview" of geoserver and it works fine with OpenLayers2.
What can it be the issue?


texts:
I have numpy array that I know bbox for and can transform like this:

import rasterio.plot
from rasterio.transform import from_bounds

t = from_bounds(minx, miny, maxx, maxy, sorted_np[3,:,:,7].shape[0], sorted_np[3,:,:,7].shape[1])

rasterio.plot.show(sorted_np[3,:,:,7],transform=t)

I want to open() this into rasterio transformed as above.

Using the comment below, tried writing new dataset as:

new = rasterio.open('new.tif',
'w', driver='GTiff',
height=sorted_np[3,:,:,7].shape[0],
width=sorted_np[3,:,:,7].shape[1],
dtype=sorted_np[3,1,1,7].dtype,
count=1, crs='+proj=latlong', transform=t)

new.write(sorted_np[3,:,:,7], 1)
new.close()

When I go:

rasterio.plot.show(new)

I get error:

TypeError: Image data cannot be converted to float

FYI:

print(sorted_np[3,1,1,7].dtype)

float32

print(sorted_np[3,:,:,7])

[[0.2168 0.2067 0.1924 ... 0.2053 0.196  0.1746]
 [0.2237 0.2218 0.2158 ... 0.2035 0.197  0.1734]
 [0.2251 0.2266 0.2251 ... 0.2018 0.1956 0.181 ]
 ...
 [0.1611 0.1637 0.1227 ... 0.1167 0.1206 0.1187]
 [0.1486 0.1661 0.1788 ... 0.1136 0.1175 0.1193]
 [0.1652 0.166  0.1952 ... 0.1141 0.1172 0.1203]]

It would be nice to not have to write a temporary tif file to do this because that seems inefficient.


texts: Querying lines that meet in endpoints to show result in clockwise order?

texts:
I have some polylines - roads. They meet in an intersection like on the picture:

I want to find out which roads meets Road A - endpoint to endpoint - and i want to know the clockwise order of the result. Thus, the example in the picture should return in the order:

{ROAD C, ROAD D, ROAD B}

Using IFeatureClass.Search() with a ISpatialFilter only returns the features ordered by OBJECTID (The road names illustrate the order og OBJECTID).

Do I have to do the math myself, or is there some smart esri woodoo that I can do?


texts: How to group Location coordinates based on Zipcode?

texts:
I'm new to Zipcode, Location Coordinates subject. I am working on a Mobile App where we are going to implement Openweather JSON.

Each Zipcode has a fixed coordinate point mentioned in the JSON.

As per my understanding, coordinates will vary on every move.

How can i group the coordinates based on Zipcode?