texts: Check geometries distorting polygons when removing gaps

texts:
Using QGIS 3.10.9 with the check geometries plugin (v0.1), and I have a simple polygon layer as a test run with a deliberate gap.

I use the check geometries plugin with the topology check 'check for gaps smaller than..' enabled and set to 100 map units sqr.

The deliberate gap is identified by the plugin but when I try to resolve it using the default method 'add gap to neighbouring polygon with longest shared edge' it distorts the shape of the polygon rather than filling in the gap.

I have previously used this plugin to fill in gaps and didn't come across this issue then.
Why is this happening?


texts: I need to label only the latest feature in a particular feature class. I actually would like to automate it, so that when I will add any new feature to this feature class, the label will automatically be assigned to the latest one. Could anyone please help me with line of code?


texts:
I am working on a story map and for my layer I would like to have a picture for every single location to be shown on the pop-up. I have followed everything on this link: http://doc.arcgis.com/en/arcgis-online/create-maps/configure-pop-ups.htm , but other than getting a picture icon in my pop-up, the picture doesn't show up.

I wonder if it is something I should change in my attribute table, like using an actual picture in my computer, or a picture on the web? or some data formatting I need to do in order to make this happen.


texts: Dissolve multiple polygons into one based on background layer

texts: I have a layer of regions (25 regions) and a layer with thousands of polygons. I want to dissolve this into one polygon for each region (so I can calculate the exact geometry). I used the dissolve tool of ArcGIS, but using this dissolves it into one layer. Is there a trick for this? Can't seem to see it in the dissolve tool itself.


texts:
I'm working in ArcMap 10.3, and I have a script that searches through all the folder and files in a directory using arcpy.da.Walk, grabs any file containing a key word, and copies it over to new a gdb based on whether it is a point, line, or polygon.

In this case I'm only grabbing vector files, either shapefiles or feature classes.

The problem is that it gets held up when it encounters non-spatial files, like a .jpg, because I'm using desc.shapeType on each file. What I'd like to know is, is there a way to quickly determine if a file is spatial or not, besides adding in some statements about its file extension (i.e.if fc.endswith(".jpg"):?


texts: Resample raster from one resolution to another with GDAL Python API

texts:
i have some trouble warping one raster to different resolution. My input is a one channel tiff with. I know the geo-coordinates of x_min and y_max. I also know that the pixel resolution is 5 meters. I now want to warp the same raster to pixel resolution of 6.5 meters.

Input in the function is the numpy array and the upper left corner coordinates with the src and dst pixel resolution.

My code is:

def resample_raster(raster, raster_x0, raster_ymax, src_xres, src_yres, dst_xres, dst_yres):

    dst_height=int((dst_yres/src_yres)*raster.shape[0])
    dst_width=int((dst_xres/src_xres)*raster.shape[1])

    dstAR=np.empty((dst_height, dst_width))

    srcDS=dstDS=gdal_array.OpenArray(raster)
    srcDS.SetGeoTransform((raster_x0,src_xres,0, raster_ymax, 0, -src_yres ) )
    dstDS=gdal_array.OpenArray(dstAR)
    srcDS.SetGeoTransform((raster_x0,dst_xres,0, raster_ymax, 0, -dst_yres ) )

    gdal_warp_options=osgeo.gdal.WarpOptions(xRes=dst_xres,yRes=dst_yres, resampleAlg=osgeo.gdal.GRA_NearestNeighbour)

    osgeo.gdal.Warp(dstDS, srcDS,options=gdal_warp_options)

    new_raster=np.array(dstDS.GetRasterBand(1).ReadAsArray())

    return new_raster

The function gives me an empty array or is killed.


texts:
I have an Oracle Spatial dataset with z-coords. I can edit the dataset in ArcMap however the z-values are always rounded to the nearest integer.

Why does ArcMap round the z-values?


texts:
I want to run the following Model

The problem is that when it reaches the copy features tool the model stops and throws an error saying that %Name%_layer (3) is not a feature class and I don't understand why.

Also if I first extract the statistics table and then combine the table with the layer in 2 different models it works fine.

Any ideas why this is happening?


texts: Can I make tiles of a raster in postgis?

texts: I am doing some work on dtms in postgis. To import the dtm, I snipped it into tiles and the work directly on the dtm works fine. Then I calculated aspect and slope of the dtm and got some bad artifacts along the tile borders, so to have this done correctly, I had to use st_union to put the tiles together before doing the calculations. So far, so good, but then I have large raster maps of slope and aspect, and when I then overlay points on those rasters, things goes very slowly (I have a point data set of a few million points). When I did the overlay on the original tiled dtm, it took about 20 minutes, when I do the same on the combined raster, it took at least one hour. Is it possible somehow to cut the calculated rasters into tiles (without exporting and re-importing them) or can I add some kind of index on those rasters that makes the overlay run faster?


texts: Deleting values from PostgreSQL database in FME?

texts:
I have a PostgreSQL database which I initialized through FME. What I want to do now is to delete some values from the database.
How can I do that? I am trying to use SQLExecutor but I don't know what to put in the input parameter called initiator.
Actually I don't quite understand why I need an input parameter. I just want to connect to the database and delete some values from the table.

What would be the workflow for this?


texts:
I am trying to let a user cancel the ArcToolbox script tool that I am creating. When the code is running, if I press cancel, the process just continues. I have looked at ESRI help on the issue (http://pro.arcgis.com/en/pro-app/arcpy/geoprocessing_and_python/understanding-cancellation-behavior-in-script-tools.htm) but I cannot get the tool to respond to the cancel command mid run.

I have tried both autoCancelling = true and false with various states of isCancelled both in and out of if statements.

What am I missing?

ArcpMap 10.3 python2.7

This is my attempt to date:

import arcpy, time
from arcpy import env

arcpy.env.autoCancelling = False

    if x == 1:
        arcpy.AddError('A landbase with this name already exists, you have 10s to cancel otherwise it will be overwritted')
        timing = 1
        while timing < 11:
            time.sleep(1)
            arcpy.AddMessage(str(timing))
            timing += 1
            if arcpy.env.isCancelled == True:
                sys.exit()


texts: Renaming fields in new layer using MakeFeatureLayer?

texts:
I have feature class that contains many fields.

I'm using MakeFeatureLayer tool in arcpy to make a new layer that would contain selected fields with different names.

Here's the code I've written to make the new layer:

arcpy.MakeFeatureLayer_management("Firebreak","Firebreaks","","","OBJECTID OBJECTID VISIBLE NONE;FIREBREAKSUBTYPE TYPE VISIBLE NONE;SHAPE Shape VISIBLE NONE;FNUMBER FNUMBER HIDDEN NONE;GPSSATAVG GPSSATAVG HIDDEN NONE;GPSCOLLECTUSED GPSCOLLECTUSED HIDDEN NONE;SURVEYER SURVEYER HIDDEN NONE;GLOBALID GLOBALID HIDDEN NONE;LVMDISTRICTCODE LVMDISTRICTCODE HIDDEN NONE;BLOCKKEY BLOCKKEY HIDDEN NONE;LASTCHANGEAPP LASTCHANGEAPP HIDDEN NONE;LASTCHANGEDATE LASTCHANGEDATE HIDDEN NONE;SHAPE.LEN SHAPE_Leng VISIBLE NONE")

The new layer contains the fields I selected, but unfortunately, those fields aren't renamed, i.e., fields of new layer have the same names as fields in original layer.

How can I make MakeFeatureLayer to rename the fields of new layer?


texts:
I want draw a corridor indicating a right-of-way around a road.
The right-of-way is 300 ft. i.e. 150 on either side of the centerline of the road.
Here I show just lines that I have hand drawn, very poorly.
How can this best be done?


texts: OpenLayers 3 modify interaction and featureover

texts:
What is the best way to add a 'featureover' event on top of an ol.interaction.Modify ?
At the map or layer level, this event is masked by the interaction.

The 'modifystart' event is fired when the feature is actually modified, that is after it is selected, and I'd like to know when the mouse is over the feature.


texts:
I am creating a model in ArcMap 10 and I need to use the Buffer Wizard tool to create multiple buffers within a polygon feature class. The problem is, I can't seem to place the Buffer Wizard within my model. I can place the Buffer or Multiple Ring Buffer command into my model, but they don't automatically populate a FromBufferDistance or ToBufferDistance value for each ring like the Buffer Wizard does. I need these values for a field calculator equation later in the model. I may be able to get by with just using the Multiple Ring Buffer tool, but it would be much more efficient to use the Buffer Wizard.

If I am unable to include the buffer wizard in my model, is there a way to run it in a script? This seems unlikely, but worth a try.


texts: ArcGIS Pro Network Analyst prioritized road map

texts:
I'm using ArcGIS Pro Network Analyst to calculate the shortest route from all points in data set A to all points in data set B, according to a road map (shp format). The road map contains two types of roads (type 1 and type 2), both in the same file but differentiated with a variable in that data called "type".

Network Analyst calculates the shortest routes using both type 1 and 2 roads. I am wondering if there is any way to prioritize type 1 roads, such that the distance is only calculated along type 1 roads unless a point is not located along a type 1 road. In that case, type 2 roads can be used to get us to the type 1-road map.


texts: Identifying layers under cursor in QGIS?

texts:
In QGIS 3.8, I've loaded hundreds of layers (georeferenced individual images), which are named by filename. I'm looking for a way to click on the map and query what image (i.e. what layer in the "Layers" panel) I've clicked on.

Is there a way to query a point in the map and have all layers underneath the point reported. There's plenty of tools for querying features, but I can't find anything for layers.


texts: Unexpected token < in JSON at position 0


texts: Where can I find a shapefile of all the Montreal Neighborhoods?

texts:
I am trying to map out Montreal and was wondering if there is already an existing polygon shapefile of all the neighborhoods.  Here is a link to a perfect example of what I am looking for.

http://en.wikipedia.org/w/index.php?title=File:ArrondissementsQuartiersMontreal.svg&page=1

If I can get that in shapefile format, it would save me a lot of time!

Thank you!


texts: How can I find the center of the circle or ellipse I just created

texts:
I need to draw circles and ellipses with a certain length/width/ radius.
I used the rectangles ovals digitizing plugin to create them, but only measured roughly the circles' radiuses. It would be great if somebody knew how to do this more exact.
what would also be good to know: is there a possibility, if I create a bigger and a smaller circle, to let the smaller circle snap to the bigger circle's center?

But the most important question is:
if I already have the circles I need, how can I make their exact center visible, so I can create my ellipses from exactly the same starting point?

I hope for help!!


texts: Raster and Geometry ST_Intersects slows when add buffer in PostGIS

texts:
I have a table containing rasters for the whole of the UK, the below query is very fast (0.1 seconds) and gets all of the rasters that are intersected by the line between a start point and an end point:

SELECT ST_UNION(rast) INTO elevationdata
FROM uk_rasters e
WHERE ST_Intersects(rast, ST_Makeline(CAST(startloc AS GEOMETRY), CAST(endloc AS GEOMETRY)));

I want to do the same, but with a buffer around the line, like so:

SELECT ST_UNION(rast) INTO elevationdata
FROM uk_rasters e
WHERE ST_Intersects(rast, ST_Buffer(ST_Makeline(CAST(startloc AS GEOMETRY), CAST(endloc AS GEOMETRY)), 75));

This was still running after 4 minutes. I assume the buffer is somehow stopping ST_Intersects from using the Index on the table. I saw mentions of ST_DWithin, but I don't see a version that takes a raster and a geometry.

Is there another approach I should be using here?


texts: For the year 2010, I want to find the percent forest coverage of several sites with a radii of 10 km around each site

texts:
Currently this is my code for generating circles on a map in GEE:

//Define all coordinates
var pointList = [
  ee.Feature(ee.Geometry.Point(-82.964662, 8.78754421)),
  ee.Feature(ee.Geometry.Point(-82.833405, 8.9453785)),
  ee.Feature(ee.Geometry.Point(-83.521707, 8.694056)),
  ee.Feature(ee.Geometry.Point(-83.33052, 8.697707)),
  ee.Feature(ee.Geometry.Point(-83.203795, 8.702278)),
  ee.Feature(ee.Geometry.Point(-83.662379, 8.685824)),
  ee.Feature(ee.Geometry.Point(-84.376866, 9.818762)),
  ee.Feature(ee.Geometry.Point(-83.485502, 8.714672)),
  ee.Feature(ee.Geometry.Point(-83.672295, 8.691173)),
  ee.Feature(ee.Geometry.Point(-83.588293, 8.694951)),
  ee.Feature(ee.Geometry.Point(-83.258253, 8.656794)),
  ee.Feature(ee.Geometry.Point(-83.06069, 8.603446)),
  ee.Feature(ee.Geometry.Point(-83.19873, 8.67821)),
  ee.Feature(ee.Geometry.Point(-82.993605, 8.785756)),
  ee.Feature(ee.Geometry.Point(-83.666655, 8.689647)),
  ee.Feature(ee.Geometry.Point(-82.976395, 8.806548)),
  ee.Feature(ee.Geometry.Point(-83.600086, 8.697011)),
  ee.Feature(ee.Geometry.Point(-83.519634, 8.562216)),
  ee.Feature(ee.Geometry.Point(-83.316672, 8.746735)),
  ee.Feature(ee.Geometry.Point(-83.518933, 8.702871)),
  ee.Feature(ee.Geometry.Point(-83.270145, 9.769966)),
  ee.Feature(ee.Geometry.Point(-83.72393, 9.16715))
];

var pointCollection = ee.FeatureCollection(pointList);

// Define a function of buffering to create circles around sites
var buffer_circles = function(feature) {
    var intermediate_buffer = feature.buffer(10000);  // buffer radius in meters
       return(intermediate_buffer); // Returns the circles
      };

// Apply function to map
var bounding_circles = pointCollection.map(buffer_circles);
    Map.centerObject(pointCollection); // Center map on sample points
    Map.addLayer(bounding_circles, {color: '#ff4c4c'}, "Buffering Circles");

I want to find the percent forest coverage of each circle for 2010. I've already looked through the GEE website, forums on this site and whatever popped up in my Googling sprees, and Youtube videos. This is my current code that I reverse-engineered using what I found for calculating forest coverage and I'm getting an error that says:

  Layer 2: Layer error: Image.projection, argument 'image': Invalid type.
  Expected type: Image.
  Actual type: ImageCollection.

Here is the code for that:

var cc = ee.Number(70)
var pixels = ee.Number(6)
var lossPixels = ee.Number(6)

//gfc2010 : global forest change in 2010
var gfc2010 = ee.Image(ee.ImageCollection("GLCF/GLS_TCC"))

var canopyCover = gfc2010.select(['tree_canopy_cover'])

var canopyCover10 = canopyCover.gte(cc).selfMask();

var contArea = canopyCover10.connectedPixelCount();
var minArea = contArea.gte(pixels).selfMask();

var prj = gfc2010.projection();
var scale = prj.nominalScale();
Map.addLayer(minArea.reproject(prj.atScale(scale)), {
    palette: ['#899fed']
});

var forestArea = minArea.multiply(ee.Image.pixelArea()).divide(10000);
var forestSize = forestArea.reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: pointCollection.geometry(),
    scale: 30,
    maxPixels: 1e13
});

I am very new to GEE and coding in Javascript.


texts: GetFeature throws IndexOutOfBounds exception

texts:
I can run a query like this using WFS 1.0.0 and 1.1.0 (the response is slightly different). But, why do I get IndexOutOfBounds exception when I try with WFS 2.0.0?

http://localhost:8080/geoserver/tiger/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=tiger:poi&maxFeatures=2&filter=<PropertyIsEqualTo><PropertyName>NAME</PropertyName><Literal>lox</Literal></PropertyIsEqualTo>


texts: Choose a Sentinel 2 image according to the amount of clouds in a given polygon - Google Earth Engine

texts:
I know that it is possible to filter an image according to the amount of clouds in the whole picture.

Is there a way to choose the image according to the amount of cloud for a given polygon instead of the whole image?

// function to create mosaic of AOI
// Mosaic the visualization layers and display (or export).
// This function clips images to the ROI feature collection
var clipToCol = function(image){
  return image.clip(geometry);
};

// Map the function over one year of data
// Load Sentinel-2
var s2 = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2017-06-01', '2017-12-30')
                  .filterBounds(geometry)
                  .map(clipToCol)
                  .sort('CLOUD_COVER')
                  .first();

//Parametros de visualización
  var rgbVis = {
  min: 0.0,
  max: 3000,
  bands: ['B4', 'B3', 'B2'],
  };

  Map.addLayer(s2 , rgbVis, 'RGB')


texts: gdal_translate to get full vector extent of geo spatial referenced PDF

texts:
I am creating GeoTif images from a USFS geospatial referenced PDFs using gdal_translate. I turned off all of the layers that I did not want and it produces a GeoTif of just the center window of the PDF. The graticule layer seems to define this.

Using QGIS, I load the PDF as a vector layer (all layers are not listed which is another question). QGIS clearly draws features that extend beyond the graticule layer. This is the PDF.

http://data.fs.usda.gov/geodata/rastergateway/data/36081/fstopo/362208152_Doe_FSTopo.pdf

Is there a way that gdal_translate can produce a GeoTif of this extra information that lies outside of the graticule layer?


texts: SQL Query error: ambiguous column name geometry

texts:
When I try a simple query

SELECT "Water_Service_Area_Boundaries_Non_Cadastral"."Network",  count(*), sum(st_length(geometry))
FROM "W_Mains_DSC_ExclAbandoned", "Water_Service_Area_Boundaries_Non_Cadastral"
WHERE ST_Intersects("Water_Service_Area_Boundaries_Non_Cadastral", "W_Mains_DSC_ExclAbandoned")
GROUP BY "Water_Service_Area_Boundaries_Non_Cadastral"."Network"

I get this error

Query preparation error on PRAGMA table_info(_tview): ambiguous column name: geometry

When I try it without the st_length it works to create a simple group with counts.

Also just the following works fine

 Select sum(st_length(geometry)), count(*)
 From  W_Mains_DSC_ExclAbandoned

Examples of the subset of the query that work independently but not when strung together

Just the total length

Just the group

When I simply add the st_length

I have tried with .geometry in the ST_Intersects and this gives correct counts but gives the pragma geometry error when the sum is added


texts: What is the modern equivalent to ArcInfo CLEAN?

texts:
What is the arcgis equivalent to the Arcinfo Workstation CLEAN command? Or perhaps better rephrased as, how does one clean topology without installing arcinfo workstation? Assume an arcinfo license level is present, just not the Arc: command line environment.

Update: A key part of the question is "without installing arcinfo workstation". I blame the confusion on ESRI's decision to use  the word "ArcInfo" for two distinct things, and then decouple them. In common usage ArcInfo can refer to  a) arcinfo workstation, and b) arcinfo license level. You can't have workstation without an arcinfo license, but you can, and increasingly more and more people do, have an arcinfo license without workstation present -- which means the answers suggesting geoprocessing clean tool can't be used.

And Kirk is right, this question can also be thought of as "What is the modern equivalent to a polygon coverage (and how do you clean it)?"


texts: Splitting lines into basic segments at vertices with PostGIS?

texts:
I have a MultiLineString table in PostgreSQL and I am trying to obtain a table with all the segments obtained from those lines if we split them at their vertices.
I have found this simple example, which actually works:

SELECT ST_AsText(ST_Split(mline, pt)) As wktcut
        FROM (SELECT
    ST_GeomFromText('MULTILINESTRING((10 10, 190 190), (15 15, 30 30, 100 90))') As mline,
    ST_Point(30,30) As pt) As foo;

But I am struggling to adapt this example to my case.

I am actually trying with something like this:

SELECT ST_AsText(ST_Split(mline, pt)) As wktcut
FROM (SELECT
     geom As mline from public.flood_def,
     (SELECT (ST_DumpPoints(public.flood_def.geom)).geom As pt FROM public.flood_def)) As foo;

where 'public.flood_def' is the table that contains the multilines. This, unfortunately, doesn't work. It throws a syntax error:

ERROR:  subquery in FROM must have an alias
LINE 30:      (SELECT (ST_DumpPoints(public.flood_def.geom)).geom As ...

Also, I'm finding it hard to understand the documentation about any of the PostGIS ST_xxxx instructions.


texts: Filtering Sentinel-5P from metadata on Google Earth Engine

texts:
I need a S5P image from a spesific date. Is there any way to filter image collection based on the metadata on GEE so the result is only from an image. I've tried something like:
var collection = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_SO2')
.select('SO2_column_number_density')
.filter(ee.Filter.eq('time_reference_julian_day', 2459022.5 ))

but it doesn't seem to work.
Which part do I need to fix?


texts: LiDAR download - ASC file to QGIS - Rasterize ->Extract -> COntour - Export to Google Earth Pro

texts:
I am trying to create a 50cm contour map using QGIS and then export it as KML to Google Earth Pro. When I get to Earth and import the KML file I am left with random red lines around the globe.
The process I am following

https://environment.data.gov.uk/DefraDataDownload/?Mode=survey

Polygon select area

Download tiles

I am selecting LiDAR Composite DTM, year 2019, resolution DTM 1M, tile ref NY96ne

That gives me files which are in various other formats but I want .asc files.
So

I have then selected LIDAR Tiles DTM 2012 1M NY96ne. Turns out they are a a series of TIFs and XML documents. Dont like the look of those.

So, I'll try LiDAR Composite DTM 2017 1M NY96ne. These show up in Explorer as Adobe Reader files, tho when opening in Adobe Reader it says the file cannot be opened as .asc is not compatible with Adobe Reader. So it appears they are .asc files

I have then opened QGIS Desktop 3.16.3 (NOT with the GRASS extension) and added in some tiles ny9568_DTM_1M and ny9569_DTM_1M, ny9468 and ny9469 and ny9668 and ny9669.

I did add a few more tiles but there are some data gaps in the data so I've removed those layers. apart from ny9469.

I have then selected all the tiles in the Layers panel and gone to Raster -> Miscellaneous -> Build Virtual Raster. This brings up the Build Virtual Raster dialogue

Up pops the Select Layers dialogue so I have selected all the layers I require

I've then set the resolution to high. This is different from my previous attempts where the res was set at average.

I then hit run, giving me a 'Virtual' layer.

14.I have unselected the other layers now, and only have the Virtual layer selected. Then I have gone Raster->Extraction->Contour and have changed the contour interval to 0.5 Hit Run

I've then right clicked on the 'Contour' layer  and then Export -> Save Feature As bringing up this dialogue box . I have added a file name and
changed the CRS from Invalid Projection to Project CPS: EPSG:4326 - WGS 84  and saved it to my desktop.

I have then opened Google Earth Pro (desktop), clicked on File -> Open -> and then selected the file I have just saved (Contours 22jan morning attempt). It opened up and the image I get is this

I accept defeat and go for caffeine. I must have spent about 20 hours on this now and I seem to get no further forward. There must be a glitch somewhere in amongst all this but that goes well beyond my capabilities on these tippytappy machines. I feel so close yet so far.

Can anyone point me in a better direction?


texts: What are the implications of invalid geometries

texts:
I've imported some data in a Postgis database and some of the geometries are reported invalid (ST_IsValidReason reports self-intersection or ring self-intersection).

The queries I am performing don't seem affected by the invalid aspect of these geometries (i'm only using ST_Distance queries).

What are the things that break when geometries are invalid?

Is fixing these geometries "automatically" (buffer(geom, 0) or ST_SimplifyPreserveTopology(geom, 0.0001)) an option?


texts: What technology to use for an Interactive web map?

texts:
Possible Duplicate:
  Comparison of JavaScript mapping libraries?

I need to build a custom interactive web map of my country with all it's regions (Peru). Which tool or technology would be the best way to go?

The idea is to have a map that shows the population of each region with a color scale: the more people the region has, the darker it's painted. Also I need to display some information as text and pictures when I click each region as the maps zooms in.

I've seen around and there are lots of "pre-developed" tools that can be customized but I need to build my own. Flash is already discarded. I was thinking of Java applet but is just an idea. Thank you for your answers.

System requirements: Apache+PHP+MySQL


texts: Map based on highest percentage from different variables


texts: CGIAR DEM converted to XYZ consists of Nodata only


texts:
I have added a layer with tracks and track points from OSM in QQGIS and exported this GPX file into an shapefile. Now I wanted to calculated the area of the buildings. I tried to do this with the tracks which also distribute the outlines of the different buidlings. But the function $area did not work in the attribute table.
So how can I get the GPX file or shape file with tracks and trackpoints into a layer with polygons?

This a picture of the track and trackpoint layer. I need the areas of all buildings to do a furter calculations.


texts: Why does ArcGIS Upgrade Geodatabase tool throw ERROR 999999: Error executing function?

texts:
I am using  ArcGIS 10.1 for Desktop to upgrade File Geodatabases using a Python script that runs fine until the step where it runs the Upgrade Geodatabase tool. The error is:

  ERROR 999999: Error executing function.\nFailed to execute

The error can also be seen when the Upgrade Geodatabase tool is run manually so I have not re-included the code which can be found in an earlier question which led to this one.

Does anyone know what would cause this error?


texts: Removing duplicate rows from a shapefile

texts:
I would like to write some c# code to delete duplicate rows from a shapefile stored in a folder.

Can anyone provide me with an example or pointers on the best way to attempt this, please?


texts: ArcCatalog File Navigation slower than Windows Explorer

texts:
I am currently using ArcMap 10.6 and have been having issues navigating to files within ArcCatalog. It has been taking me about 10-15 seconds to open folders in ArcCatalog when it only takes a fraction of a second to open the same folder in Windows Explorer.

Does anyone know the reason for this lag?


texts: Describing/creating VRT/CSVT files of variable length space separator for gdal_rasterize command

texts:
I have LAS points (DTM) files with space separators. And the problem is for Z coordinate:

 360020.678  6245020.453   0.008
...
 364895.678  6245245.453  14.970
...
 364370.678  6246870.453 112.312

String looks like: [space]X[space x2]Y[space xn]Z
Format for Z could be [space][space][space]Real(5.3), [space][space]Real(6.3) or [space]Real(7.3).
I tried with *.csvt file, but it didn't work out.

"Real(11.3)","Real(13.3)","Real(8.3)"

ogrinfo output:

field_1: Real (11.3)
field_2: Real (13.3)
field_3: Real (8.3)
field_4: String (0.0)
field_5: String (0.0)
field_6: String (0.0)
field_7: String (0.0)
OGRFeature(T_test):1
  field_2 (Real) = 360020.678
  field_4 (String) = 6245020.453
  field_5 (String) =
  field_6 (String) =
  field_7 (String) = 0.008
  POINT Z (360020.678 6245020.453 0)
...
OGRFeature(T_test):2095
  field_2 (Real) = 362370.678
  field_4 (String) = 6245270.453
  field_5 (String) =
  field_6 (String) = 15.712
  POINT Z (362370.678 6245270.453 0)
...
OGRFeature(T_test):17110
  field_2 (Real) = 362745.678
  field_4 (String) = 6247145.453
  field_5 (String) = 128.325
  POINT Z (362745.678 6247145.453 128.325)

GDAL version: GDAL 2.4.0, released 2018/12/14
OS Win10.

Is it possible to rasterize such points with the standard bin version?


texts: Geoserver API - Move layer to different workspace

texts:
Is it possible to move raster layer with coverage store to different workspace using geoserver rest api?

I suspect that i can remove layer and publish it again but maybe there is a simpler way?


texts: ESRI Web AppBuilder - custom attribute - only show fields that contain values

texts:
Is there a way to only show fields that have a value, i.e. are not NULL, in the Web AppBuilder Web Map pop-up? If this can't be done in the Web Map can it be done in the mxd attribute table?

I have a service that contains a layer with many fields but some features in the layer may not have data in all fields. To increase the visual and usability of the viewer I would like to only show fields that have a value for each feature identified on.


texts: Simplest way to display a new marker on a map already displayed?

texts:
Could someone provide an example of the most dirt simple method to force a new marker to pop up based on a given lat/lng?

I do not need js to write to my db, as that will be done server side; but I do want to display the new marker without refreshing the whole page.  The UI will give a text input that gets geocoded and inserted into the table (e.g., string that will be located on google maps "123 seaside, any town, arzonia " . I just need that new point to be marked on the map as it stands.  Again that client-side point won't be written anywhere, it's just added temporarily (until the map is hard refreshed again).

see http://highrank.robh71.com/map1/index.php

This page primitively demonstrates what I want to accomplish.  Enter a new point or location, marker shows up without a reload.  Geolocation already works, this is simply a display question.  Currently I process data control via the sql api on a separate page for simplicity.  This will obviously be changed!  I'm really only concerned with the marker rather than the data update so let's pretend that the data write doesn't even matter.).


texts:
I'm working on a modeler project. I'd like to automate the creation of web maps with qgis2web, but I can't get the algorithm 'Export vector layer' or 'Export project' work in the processing toolbox.
Have anyone tried / succeeded?


texts: How to write features in file geodatabases (FileGDB) with GDAL and Python?

texts:
I am trying to copy features from a spatiallite database to a filegeodatabase. The target FileGDB already exists with the required feature classes.
I am using the FileGDB-API 1.4 on Ubuntu and Python 2.7.6.

from osgeo import ogr
#define db paths
psimgdb='/home/ubuntu/Simulation.gdb'
psrcdb='/home/tmp/2011.db'

#open connections
indriver=ogr.GetDriverByName('SQLite')
srcdb = indriver.Open(psrcdb,0)
outdriver=ogr.GetDriverByName('FileGDB')
simgdb = outdriver.Open(psimgdb,0)

#get input and output layers
inp_main=src.GetLayer('conduit')
out_main=gdb.GetLayer('msm_link')

#cycle through the features in the input layer
for feature in inp_main:
    outfeature = ogr.Feature(out_main.GetLayerDefn())
    outfeature.SetGeometry(feature.GetGeometryRef())
    outfeature.SetField('ID',feature.GetField('id'))
    out_main.CreateFeature(outfeature)
    outfeature.Destroy()

#close connection
simgdb=None
srcdb=None

The script runs fine, when I debug the code I can see that the attributes are created correctly on the target features, but no data is written to the file geodatabase "simgdb". Am I missing something here?
I have also tried del out_main, simgdb.Destroy(), del simgdb and

out_main.StartTransaction()
out_main.CreateFeature(outfeature)
out_main.CommitTransaction()

as shown for PostGIS databases in the Python GDAL/OGR cookbook, without success.

--UPDATE--
I actually realized that the CreateFeature function returns with error code 6. Not sure what this means...


texts: How to work with OS map files from Digimap?

texts: I am trying to add map tiles as raster files in QGIS but when I add them they all appear on top of each other. I did this once before and everything was ok but I am now unable to find how to get the tiles showing as one map. I am using OS Streetview.


texts: I am a beginner in QGIS. I am trying to create 3D models using Qgis2threejs plugin. I installed it but am not able to find the tool. I only see "Qgis2threejs Exporter/Exporter without Preview tool" under Web.


texts: Is there a documented convention for the longitude of the North and South poles?

texts:
I am working on plus codes, which is a formal specification for encoding areas bounded by north, south, east and west latitudes and longitudes.
We need to decide if the North and South poles shall be specified as belonging to one or infinitely many possible encodings.
The mathematician may say "they're just numbers, encode (+90, X) however you want". But the application designer might say "every place should encode to one coordinate and every coordinate to one place".
Is there any geographic information system specification which assigns the North and South 90th "parallel" to a specific longitude?


texts:
I'm using QGIS to extract coordinates of a road network.
Concretely, I downloaded a mbtiles file from Internet and loaded the file into QGIS like following:

Then I tried to export the layer transportation to a csv file. However, I didn't see the coordinates (latitude, longitude) of the road.

I tried to add the coordinates by using Vector > Geometry Tool > Export\Add geometry column, but it doesn't produce the coordinates.

Furthermore, in the attribute table, I can see _row and _column.
How can I obtain the latitudes and longitudes for the road system?


texts: DEFINE_SHAPE Global Mapper Clip/Crop Raster on Import


texts: Limit Geocode results based on application extent?

texts:
I am creating an Esri Flex API application (not the viewer), and have been creating a tool to use the World Places Locator.

For now, I am just zooming to the top scoring candidate returned.

If you search for Brighton, it returns you 20 candidates, sorted by score.

I read on an older blog post:

  You can also filter the results by extent using client logic. The
  candidates field to look for in this case are North_Lat, South_Lat,
  East_Lon, and West_Lon.  Source.

I am thinking it would be a good idea to have a tickbox on my search, that limits the results to the current extent of the map in the Flex app.  (Or maybe buffer a litte out from the view).

I understand the logic on how to go about this, but was wondering if anyone has already done this (does not have to be Flex) so that I can save some time coding it.

If not, I will have a go at this next weekend and post my answer.


texts: Sort attribute table in order of selection

texts:
I have a large number of planned borehole sites. I want to label them in the Attribute Table in order of priority. I would like to select them in that order (priority) and have the Attribute table list them in the order of selection so that I can quickly enter text order.

I can't seem to do this as Show Selected Features doesn't list them in the order selected?

I am using QGIS 3.4.


texts: Importing matplotlib basemap and shapely

texts:
I'm not sure if this is the right place to ask this, but here we go.

I have a python script that uses matplotlib's basemap and another part that uses shapely to do an intersection of 2 polygons.  If basemap is imported before shapely and I run the intersection I get this exception:

  intersect_poly = grid_poly.intersection(data_poly)
File "/sw/lib/python2.7/site-packages/shapely/geometry/base.py", line 334, in intersection
  return geom_factory(self.impl['intersection'](self, other))
File "/sw/lib/python2.7/site-packages/shapely/topology.py", line 53, in __call__
  "This operation produced a null geometry. Reason: unknown")
shapely.geos.TopologicalError: This operation produced a null geometry. Reason: unknown

If I import shapely first everything works fine.  I would assume this is because of some "funkiness" in the way they are accessing the GEOS library.  I've checked that in both situations the same library file is loaded in shapely ("print shapely.geos._lgeos").

Does anyone have an idea as to why this is happening and if there is a right way of doing this?  Does this happen for anyone else? In the mean time I can just make sure to import shapely first (not sure if that affects basemap yet).  Otherwise maybe I'll skim through the basemap source and ask this question on one of the mailing lists.  Any suggestions on which mailing list I should post this to? Thanks.

I'm using OSX(10.7) with a fink install that has "libgeos3.3.3-shlibs", "libgeos3.3.1-shlibs", "libgeos3.3.1", "libgeos3.3.0-shlibs", "libgeos3.3.0", and "shapely-py27 (1.2.16-1)" installed.

And here's a simple test script that reproduces the problem (flip the imports and it works):

from mpl_toolkits import basemap
from shapely import geometry

g_ring = [(-88.462425, 26.992203), (-57.847187, 26.992203), (-57.847187, 17.599869), (-88.462425, 17.599869), (-88.462425, 26.992203)]
grid_g_ring = [(-123.044, 59.844000000000001), (-49.384999999999998, 57.289000000000001), (-65.090999999999994, 14.335000000000001), (-113.133, 16.369), (-123.044, 59.844000000000001)]

data_poly = geometry.Polygon(g_ring)
grid_poly = geometry.Polygon(grid_g_ring)

print grid_poly.intersection(data_poly).area


texts:
I am using postgresql 9.4.1 and postgis 2.1(pg 94) and also qgis 2.4.0. All the installation went well. But, when I tried to import shapefile to postgis via postgis shp & dbf loader it couldn't work! Here is the warning message:

  ==============================

  Importing with configuration: area, public, geom, D:\temp\area.shp,
  mode=c, dump=1, simple=0, geography=0, index=1, shape=1, srid=0
  Shapefile type: Arc PostGIS type: MULTILINESTRING[2] Failed SQL
  begins: "SET CLIENT_ENCODING TO UTF8; SET STANDARD_CONFORMING_STRINGS
  TO ON; BEGIN; CREATE TABLE "public"."area" (gid serial, "objectid"
  int4, "fid_block3" int4, "shape_leng" numeric); ALTER TABLE
  "public"."area" ADD PRIMARY KEY (gid); SELECT AddGeometryColumn('pu"
  Failed in pgui_exec(): ERROR:  fungsi addgeometrycolumn(unknown,
  unknown, unknown, unknown, unknown, integer) tidak ada LINE 9: SELECT
  AddGeometryColumn('public','area','geom','0','MULTILI...
                 ^

  HINT:  Tidak ada fungsi yang sesuai dengan nama dan tipe argumen yang
  diberikan.Anda mungkin perlu menambakan cast tipe eksplisit

  Shapefile import failed.

  ==============================

At first I thought it was because of my shapefile is broken. But it was not. I try to produce new shapefile and still couldn't import!
If anyone ever found the same case, please share your solution to me.

Oh ya, in addition, I also have tried the SPIT tools in qgis and still I can't import the shapefile with warning massage:

Problem inserting features from file:
D:/temp/for_postgres/poligon.shp

Error while executing the SQL:

  SELECT f_table_name FROM geometry_columns WHERE
  f_table_name='poligon' AND  f_table_schema='public'The database
  said:ERROR:  relasi Â«  geometry_columns Â» sudah ada LINE 1: SELECT
  f_table_name FROM geometry_columns WHERE f_table_name...^


texts: How to do Reverse Geocoding from MySQL/PostGIS Sql ???

texts:
I want to do reverse geocoding from MySQL or PostGIS SQL.

I have shapefile up to town level of India. So I have a geometry boundary for each town.

What I want is to first find the nearest point from the current position and check that point is in which towns.

I want to get the result in this format: Towns, District, State, Country

I dont have knowledge of PostGIS SQL but I know MySQL. But if anyone can help me I will use the PostGIS SQL also.

So I want fastest way to do reverse geocoding.

Thank you.


texts: Get coordinates from web mercator projection to Lan & Lot (JavaScript/ node.js)

texts:
I have the WKT coordinates like below.

LINESTRING (16139582.976129573 -4541142.959679384,16139974.715899538 -4540837.211566243)

How can I get the Longitude and Latitude from that.

I am using GISCloud.


texts:
Based on this link example, I am trying to use Point feature but I get :

my.js:10 Uncaught TypeError: ol.source.Point is not a constructor

   var iconFeature = new ol.Feature({
        geometry: new ol.source.Point([35.1495, 33.39243]),
        name: 'A point',
        population: 4000,
        rainfall: 500
    });


texts: ArcGIS ModelBuilder Iterate Tables and XY Event Layer?

texts:
I'm running into an issue that was not addressed in a previous question, but linked to a similar issue (ArcGIS modelbuilder: Iterate through Excel sheets?).

I'm using ModelBuilder to import .xlsx tables into Arc and turn them into point data.

My problem is not based in the naming of the outputs or converting the points to feature classes,  it's actually linked between the Iterate Table tool and Make XY Event Layer tool.

Whenever I connect my excel workbook to iterate tables, the output table is directly connected to the make XY event layer (as seen in the first answer on the link) but does not seem to actually connect the table coordinates. In Make XY Event it has no X or Y fields to display. I checked on Esri's page on the tool and my data meets all the requirements for the tool. I even ran the tables separately in the tool to verify that they would work. Everything runs fine until the incorporation of the Iterate Tables tool. At one point it actually connected with no problem and then disconnected when I tried to verify the model. I can't seem to find the answer to why this is happening on google either. I'm fairly new to ModelBuilder so I'm afraid I'm just making a simple mistake with the iterate tool, but there isn't much to actually mess up so I'm stumped.

Extra details, the table used only contain Lat and Long data to minimize issues and I'm using ArcMap 10.2.2


texts: Identifying upside and downside of road using QGIS?

texts:
I am a forest engineer and I have to reclassify a DTM in upside pixels and downside ones in relation to a line shapefile indicating skid trails network of the forest parcels.
I don't have necessarily to reclassify the DTM but I need to know and identify zones which have an higher altitude than the skid trails and zones which have lower one.

I am using QGIS.


texts: Deleting portion of raster using ArcMap?

texts:
How do I delete a portion of a raster image using ArcMap?

Attaching the raster image and the portion to delete.


texts: How do I join ACS data to this Census shapefile?

texts:
I'm using QGIS.  I am trying to join the .CSV data from the American Community Survey (ACS):
https://gisdata.mn.gov/dataset/us-mn-state-metc-society-census-acs
With the shapefile from the 2010 Census:
https://gisdata.mn.gov/dataset/us-mn-state-metc-society-census2010realign
Both files have GEOID in common.  But in the ACS file GEOID is a longer number that includes letters, unlike in the other file, so it doesn't appear to refer to the same thing.  It does, however, have a field called GEOID2 that appears to be the same as the GEOID in the other file. Should I change the title of GEOID2 to "GEOID" and use this as a common field to join them?  Or is there a better way to do this?


texts: How to set “draw line inside polygon” using Python?

texts:
In the QGIS properties of a layer I want to activate "draw line only inside polygon" for a specific layer, but using Python. I don't want to use the GUI for it.

Can this be done?

I set the line width with the following code:

props['line_width'] = '1.0'
lineLayer = QgsSymbolLayerV2Registry.instance().symbolLayerMetadata("SimpleLine").createSymbolLayer(props)

Does anyone have any idea how to check the "draw line only inside polygon" checkbox?

Please see attached image.


texts: Using ArcPy to Import Symbology instead of Apply Symbology?

texts:
I am using arcpy (ArcMap 10.1) to create a series of maps of the same extent that show increased count of disease in areas over time. For example, in the first map, the disease count will be low in each area, the second map will be slightly more, etc. I created a layer file with a symbology that I want to use in each map so that, as time increases, people can see how the count of disease is increasing in each area relative to prior maps.

In ArcMap itself, I can go to the Layer Properties and Import

my layer file so that ranges for each coloring bin do not change. However, it seems that arcpy only has access to the Apply Symbology From Layer tool. This applies the rules of the layer file I have (ie. Natural Breaks, 10 classes), but changes all of the bin sizes and makes them only appropriate to that specific layer.

Is there any way to access something like the Import button from ArcPy?


texts:
Is there a way to add rows to an existing dataset while importing in a table?

I've have tried inserting directly from an existing dataset but get an issue withe cartodb_id already existing.


texts: Split field name of feature class in ArcGIS

texts:
Is it possible to split a field name within an attribute table of a feature class into two new fields where just the field name is being split. Please note that I do not need to split the features (rows), just the field name.

So for example, I have a field name called "JoinCount_5meters" and I would like to split this field name into two new fields called "JoinCount" and "5meters". Again, the features themselves (row values) do not need to be split. so, anyone have any ideas how I could split a field name into two new field names in ArcGIS?

I understand I can do this manually but I have 250 feature classes for which I need to split 16 field names in each feature class. It would be great if I can do this process using Python or within model builder using an iterator.


texts: Writing raster with attribute table in R?

texts:
I fitted a species distribution model and predicted the suitable habitat using this model, and then I got a map plotted in R.

I tried to save this prediction as a raster, using code:

rf<-writeRaster(px, filename="predict.tif", format="GTiff")

However, I could only export a tif format file without an attribute table. I could put this raster to ArcGIS software, but without an attribute table, it is difficult to run any analysis.

How can I export a raster with it's attribute table in R?


texts: GeoTIFF / Geotools / Mollweide / customer projection not recognized

texts:
Although I got a hint from Ian, I still have a problem to read a GeoTIFF in Geotools.
I ask this as a new question so that I can put all info.

The GeoTIFF is from
http://cidportal.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_POP_GPW4_GLOBE_R2015A/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k/V1-0/

gdalinfo GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0.tif
Driver: GTiff/GeoTIFF
Files: GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0.tif
       GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0.tif.ovr
       GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0.tif.aux.xml
Size is 35497, 15236
Coordinate System is:
PROJCS["World_Mollweide",
    GEOGCS["GCS_WGS_1984",
        DATUM["D_WGS_1984",
            SPHEROID["WGS_1984",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Mollweide"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",0.0],
    UNIT["Meter",1.0]]

As it is a Mollweide projection, I add the EPGS definition.

epsg.properties file:

54009=PROJCS["World_Mollweide",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mollweide"],PARAMETER["Central_Meridian",0.0],UNIT["Meter",1.0]]

    URL epsg = getClass().getClassLoader().getResource("epsg.properties");

    Hints hints = new Hints(Hints.CRS_AUTHORITY_FACTORY, PropertyAuthorityFactory.class);
    ReferencingFactoryContainer referencingFactoryContainer = ReferencingFactoryContainer.instance(hints);
    PropertyAuthorityFactory factory;
    factory = new PropertyAuthorityFactory(referencingFactoryContainer, Citations.fromName("EPSG"), epsg);
    ReferencingFactoryFinder.addAuthorityFactory(factory);
        ReferencingFactoryFinder.scanForPlugins();

This prints correctly:

System.out.println(CRS.decode("EPSG:54009"));

PROJCS["World_Mollweide",
  GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",
      SPHEROID["WGS_1984", 6378137.0, 298.257223563]],
    PRIMEM["Greenwich", 0.0],
    UNIT["degree", 0.017453292519943295],
    AXIS["Longitude", EAST],
    AXIS["Latitude", NORTH]],
  PROJECTION["Mollweide"],
  PARAMETER["semi_minor", 6378137.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  PARAMETER["central_meridian", 0.0],
  UNIT["m", 1.0],
  AXIS["x", EAST],
  AXIS["y", NORTH],
  AUTHORITY["EPSG","54009"]]

But

String fileName = "/home/XXX/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0.tif";

new GeoTiffReader(new FileInputStream(new File(fileName)));

still throws an error

GeoTiffMetadata2CRSAdapter::createCoordinateSystem:Only Geographic & Projected Systems are supported.
    at org.geotools.coverage.grid.io.imageio.geotiff.GeoTiffMetadata2CRSAdapter.createCoordinateSystem(GeoTiffMetadata2CRSAdapter.java:186)
    at org.geotools.gce.geotiff.GeoTiffReader.getHRInfo(GeoTiffReader.java:329)
    ... 4 more

geotools version 21.0
maven

<dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-main</artifactId>
            <version>${geotools.version}</version>
        </dependency>

    <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-referencing</artifactId>
            <version>${geotools.version}</version>
        </dependency>

        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-epsg-wkt</artifactId>
            <version>${geotools.version}</version>
        </dependency>

    <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-geotiff</artifactId>
            <version>${geotools.version}</version>
        </dependency>


texts: DC.js - barChart not filtered on click of geoChoroplethChart

texts:
I have a choropleth map of counties, and a bar chart showing population of each county on the map. When I click on the bar chart the map is filtered and redrawn to show the selected bar, but when I click on the a county in the choropleth map the bar chart is not filtered to show the population data. I don't understand why it would be filtering one way but not the other.

    Population by counties (color: total population)
    reset
     | Current filter:

    Population by county
    reset
     | Current filter:

    Reset All

var iowaMap = dc.geoChoroplethChart("#iowa-map");
var populationChart = dc.barChart("#population-chart");

d3.csv("iowaCountiesPop.csv", function (data) {

      data.forEach(function (d) {
        d.county = d.county;
        d.popByCounty = +d.e2015;
      });

      var data = crossfilter(data);

      var counties = data.dimension(function (d) {
        return d.county;
      });

      var counties2 = data.dimension(function (d) {
        return d.county;
      });

      var popByCounty = counties.group().reduceSum(function (d) {
        return d.popByCounty;
      });

    d3.json("IowaCounties.json", function (countiesJson) {
        iowaMap.width(990)
                .height(500)
                .dimension(counties)
                .group(popByCounty)
                .projection(d3.geo.mercator()
                    .translate([495, 250])
                    .rotate([93 + 20 / 60, -41 - 60 / 60])
                    .scale(7900))
                .colors(d3.scale.quantile().range(colorScheme[quantiles]))
                .colorDomain([0, 430640])
                .overlayGeoJson(countiesJson.features, "NAME", function (d) {
                    return d.properties.NAME;
                })
                .title(function (d) {
                    return d.key + " County \nTotal Population: " + numberFormat(d.value);
                })
                .on('renderlet', function(map) {
                    map.selectAll("path").on("click", function(d) {
                    //console.log("click!", d)
                     map.filter(d.properties.NAME)
                        .redrawGroup();
                    })
                });

         populationChart.width(width)
            .height(height)
            .dimension(counties2)
            .group(popByCounty)
            .x(d3.scale.ordinal().domain(counties))
            .xUnits(dc.units.ordinal)
            .margins({top: 0, right: 0, bottom: 70, left: 70})
            .yAxisLabel(["Population Values"])//,[12])
            .xAxisLabel("County Names")
            .barPadding(0.1)
            .outerPadding(0.05)
            .elasticY(false)
            //.turnOnControls(true)
            .on('renderlet', function(chart) {
                chart.selectAll('rect').on("click", function(d) {
                    //console.log("click!", d)
                    chart.filter(d.data.key)
                        .redrawGroup();
                })
                chart.selectAll("g.x text")
                    .style("text-anchor", "end")
                    .attr("dx","-8")
                    .attr("dy", "5")
                    .attr("transform", "rotate(-50)");
            });

        dc.renderAll();

    });

});


texts: Replace geometry of polygon feature using UpdateCursor and SHAPE@ token - arcpy

texts:
Within a Python script, I would like to replace the geometry of all features in Layer A with the geometry of features in Layer B if their ID fields match.

I am not sure how to do this, but the closest I have gotten (conceptually) is by using updateCursor and the SHAPE@ token.  The following code kind of makes sense to me, but is not updating the geometry.  Does anyone have any suggestions for fixing this code or another method that would work better?

def burnCPP(sgcn_layer, cpp_layer):

    # create cursors for Layer A and Layer B
    sgcn = arcpy.da.UpdateCursor(sgcn_layer, ["DataSource", "DataID", "SHAPE@"])
    cpp = arcpy.da.SearchCursor(cpp_layer, ["EO_ID", "SHAPE@"])

    # create list and fill with IDs from Layer A
    sgcn_list = []
        if row[0] == "PNHP":
            sgcn_list.append(str(row[1]))

    # create list and fill with matching IDs and geometry of Layer B
    match = []
    if str(row[0]) in sgcn_list:
        match.append([row[0], row[1]])

    # start iteration through rows in Layer A
    with sgcn as cursor:
        for row in cursor:
            # for each row in Layer A, start iteration through each list in match
            for lists in match:
                # if ID in row in Layer A matches ID in list in match fill geometry in Layer A with geometry stored in match list
                if row[1] == lists[0]:
                    row[2] = lists[1]
                    cursor.updateRow(row)


texts: Overlaying Line (Vector) Layer on Raster Zones?

texts:
I am using ArcGIS 10

I have a raster layer with 3 values that make up zones. For the purpose of this example, lets say these zones are forest, grass, and open water. My whole raster is made up of these three values. I also have a line layer (it is currently a vector file but I'm fine converting to raster if I need to). This line layer is streams. I want to attach this stream layer to the raster layer based on zones. so I will end up with 6 values (Forest with stream, forest without stream, grass with, grass without, open water with, open water without). I basically want to do an overlay, but apply the stream value to the whole zone and not just a pixel. If a stream falls into a forest pixel, I want all neighboring forest pixels to also have that stream.

I really want to avoid converting the raster to vector. I looked into the zonal function tools but they seem to only give me statistics and values based on zones rather than overlay two layers based on zones.


texts: Docking Symbology / Layer Properties window in QGIS?

texts:
In QGIS, is there a way to make the "Layer Properties" (aka the "Vector Properties" window) docked to the right side of the screen?

More specifically, the "Symbology" tab.

I ask this because sometimes I just want to more easily change the reference attribute the map is being colored by, and having to double click on the layer to open the "Layer Properties" window and choosing the correct tab adds a few pesky clicks I would like to avoid.

I found an answer at Displaying QGIS attribute table on the right for another window, but I'm having a hard time adapting the answer to this case.

I'm currently using QGIS 3.2.


texts: How to speed up work with local maps in QGIS?

texts:
within one project we are distributing a QGIS project including offline maps (you can get it here - has about 1GB) to display measurement data. Looks this way:

I need to keep it as easy as possible so I used spatialite/sqlite DBs with custom styled OSM data (Google like styles from @Underdark) and contour lines generated from EU-DEM:

However the performance is relatively slow although I set particular layers to only show in defined scale range. I can avoid much lagging by first switching off the offline maps, then add new data, zoom to new layer and switch offline maps again on but it is not much comfortable.

I cannot use additional SW like PosGIS DB as it would be complicated for the users to install it. Is there any way to improve the performance?

I am thinking about generating raster layers in 1:200 000, 1:50 000 and 1:25 000 scale and using maybe virtual raster + scale dependent display to avoid forcing the user to manually switch the offline map layers.

What would you recommend?


texts: Call ArcGIS Server Buffer command without using SOM

texts:
I wanted to see what the best way was of calling the geometry service to carry out a Buffer operation on ArcGIS Server 10, without using the SOM.  Ideally I want to consume the REST API, using a wrapper so that it was easy to call from c#.

Is the GISClient the way forward?  From what I can determine, it looks like you only get stateless (and thus low footprint) access to the Mapserver and Geocode services.

I know that I could construct a call using a URL, but I wondered if there was C# API that wrapped the REST API for you.


texts:
I have imported a shapefile whole polygon is diagonally oriented. Rather than rotate the polygon, I would like to adjust the scale bar/legend/title to fit this unique orientation.

Would like to know it both through GUI and python.


texts: If/then python for classifying existing field values to a new field

texts: I'd like to have a tool or quick python script to classify an existing number field to a new text field.  For example I have a roads layer that has a type identifier as a coded number and I would like to create a new field that described the number code as text.  So say the existing field 'FCODE" had code "210" , and I want to write "Highway" in a new field "Type" and do that for multiple number-text pairs. I could select by attribute and classify that way but it would be nice to do it all in one step.


texts: QGIS - Orfeo Toolbox: how to run support vector machine classification?

texts:
I am trying to classify SPOT images with the support vector machine classification in Orfeo Toolbox with QGIS. I installed Orfeo Toolbox with OSGeo4W. I have made a ROI image with 600 training samples (2 classes), and used TrainImagesClassifier (svm). However when I use Image Classification to classify the image with the svm, it gives me an error. The log only says:

  the following layers were not correctly generated: Output image,
  Confidence map

I have tried using different images, leaving it to save to temporary file, and changing the save directory but nothing works.


texts: How to use LIKE clause/condition in Arcpy.Exists/Delete

texts:
I would like to delete several raster files in an ArcSDE SQL Server database.  I would like to do so with Arcpy.  How can I succesfully use a LIKE condition to delete all files in my database that have a certain string of characters, such as _.  Meaning that I want to delete all files that contain an underscore.  Using the example from Esri, where would I delete all files that for instance contain roadbuffer.

import arcpy

# Set the current workspace
#
arcpy.env.workspace = "c:/base/data.gdb"

# Check for existence of data before deleting
#
if arcpy.Exists("roadbuffer"):
    arcpy.Delete_management("roadbuffer")


texts: Converting state plane coordinates to latitude/longitude

texts:
I'm trying to convert pairs of x,y coordinates to latitude/longitude using python.  The projection is "NAD_1983_HARN_StatePlane_Oregon_North_FIPS_3601_Feet_Intl" and a sample x,y pair is 7625659, 714714.  So far, I have found this question and have tried using the code in the first solution:

from pyproj import Proj, transform

inProj = Proj(init='epsg:????')
outProj = Proj(init='epsg:????')
x1,y1 = 7625659,714714
x2,y2 = transform(inProj,outProj,x1,y1)
print x2,y2

However, I'm having trouble finding an epsg code for the location.  I came across this but I can't seem to find the epsg code for that area.  All I can see is the sr-org code.  Also, it seems strange that the location on the map is not in Oregon.  I was able to convert it (the output looked reasonable) using the tool at http://beta.ngs.noaa.gov/gtkweb/ using datum=NAD83, projection for conversion=SPC, Norting=714714, Easting=7625659, Units=International feet and SPC Zone=OR N-3601.  How should I convert this projection to latitude and longitude?

Edit:
I tried using the tool at http://prj2epsg.org/search which gave me 2913 (I also tried 3646 per BradHards suggestion, it came up the same) for inProj and 4326 (BradHards) for outProj, but it's giving me a location (-62.0472510133 32.7243167795) which is in the middle of the ocean, not in Portland.  The lat/lon from http://beta.ngs.noaa.gov/gtkweb/ is (45.605018016, -122.752592697) which is in Portland (Sorry, I can't comment on the site...).


texts: Understanding Coordinate System Transformation in ArcMap?

texts: For the Data Frame Properties on a map in ArcMap, what is the difference between selecting a new coordinate system and using a Transformation?


texts:
I am trying to mirror or rotate(180) my layer some thing like  this "-->" to this "<--".
Is there any built in plugin or support available to do this job or i have to write code in pyqgis?

I don't know what the exctly term is better to asked either it is invert, filp, rotate or mirror.


texts:
I have a base layer on top of which I will plot two separate feature layers, one containing polygons and the other containing points. Now I have to find the points which plot into the polygon. I tried various examples like  https://developers.arcgis.com/javascript/jssamples/util_relation.html but somehow things are not getting done.

Can anyone suggest a better idea?


texts: Where_Clause expression working for one FeatureClassToFeatureClass_conversion export but not another

texts:
I am attempting to export a feature class to both a shapefile and feature class with a where_clause restriction on particular values within particular fields. I have written a function to take dataset paths and export subsets for particular datasets.

I am getting correct/incorrect outputs depending on the particular dataset that is passed to the function. I have triple checked field names/values and implemented AddFieldDelimiters to try to circumvent dataset format differences.

The code can be viewed below:

def export_data_retain_fields(dataset, output_directory):
    print("Exporting {} to {} directory".format(dataset, output_directory))
    dataset_fname = DCM_sde_new_names[dataset.split('.')[-1]].replace('DCP_TRD_', '')

    '''
    Below section of code doesn't work. Defining field names in sql
    expressions require specific formatting depending on if
    the data is in shapefile or feature class format. The AddFieldDelimiters
    method allows you to not have to format your
    sql expression to account for if it's a shapefile or feature class
    '''

    if dataset == DCM_sde_path:
        print("Restricting Feature Type and Type fields for
        {}".format(dataset_fname))
        sql_expression = """Type <> 'Former_ROW' AND Feature_type <> 'Demapped_St'"""

        # Here is where I declare export with where_clause equal to sql_expression defined above. This DOESN'T work

        arcpy.FeatureClassToFeatureClass_conversion(dataset, output_directory,
                                                 dataset_fname, where_clause=sql_expression)

    '''
    Below section of code works. Defining field names in sql expressions
    require specific formatting depending on if
    the data is in shapefile or feature class format. The AddFieldDelimiters
    method allows you to not have to format your
    sql expression to account for if it's a shapefile or feature class
    '''

    if dataset == DCM_SCL_sde_path:
        print("Restricting Route Type and Feature Status for
        {}".format(dataset_fname))
        sql_expression = """Route_type <> 'Former_ROW' AND Feature_status <> 'Demapped'"""
        print(sql_expression)

        # Here is where I declare export with where_clause equal to sql_expression defined above. This DOES work

        arcpy.FeatureClassToFeatureClass_conversion(dataset, output_directory,
                                                dataset_fname, where_clause=sql_expression)
    else:
        arcpy.FeatureClassToFeatureClass_conversion(dataset, output_directory,
                                                dataset_fname)

The function works correctly for the DCM_SCL_sde dataset path but returns the full dataset (omits the subset) for the DCM_sde dataset path. I have triple checked field names and values and have not been able to find any discrepancies.

I have ran the same exact query directly in ArcMap and got a correct data subset to export. I also attempted to use Select_analysis with arcpy rather than FeatureClassToFeatureClass with no success.


texts: Renaming outputs with iterative ModelBuilder?

texts:
My model is meant to create mosaic datasets within a gdb. I am using the "For" function, and setting it to 6. The output correctly creates 6 mosaic datasets, each with the prefix T.

I would like these to have a different name, grabbing the name of the gdb and then _# (so G_057329641020_01_1, G_057329641020_01_2, etc.). How can I add the renaming function?


texts: Creating an sf object from the maps package

texts:
I am trying to use data from the maps package in an sf workflow. I'm not sure if I have the operation correct (see below). It seems to work until state 20 (massachusetts). What am I doing wrong? Maybe there is a better way to do this?

state_data <- ggplot2::map_data("state")
state_data_sf <- lapply(unique(state_data$region),
                function(x)
                  list(matrix(unlist(
                  state_data[state_data$region == x,][,c("long", "lat")]),
                  ncol = 2)))

sf::st_polygon(state_data_sf[[19]])

  POLYGON((-87.4620056152344 30.3896808624268, -87.4849319458008 39.6200332641602, -78.1113357543945 39.6658706665039, -78.1342544555664 39.67732...

sf::st_polygon(state_data_sf[[20]])

  Error in MtrxSet(x, dim, type = "POLYGON", needClosed = TRUE) :
    polygons not (all) closed


texts: Table to Table script fails when copying to SDE

texts:
Since TBL2SDE no longer exists (we're running a 64-bit server and ArcGIS 10 doesn't include that executable on a 64-bit installation) I've been using the Table to Table script in ArcToolbox to copy some Info and DBF tables to SqlServer. This has been working fine but I ran into a problem today where it simply does not allow me to copy any tables.

When I run the command, it returns:

ERROR 999999: Error executing function.
Item with the same path name already exists.
DBMS table not found [Microsoft SQL Server Native Client 10.0: Invalid object name 'gismo.RFCGIS.aocondo'.] [gismo.RFCGIS.aocondo][STATE_ID = 19653]
DBMS table not found [Microsoft SQL Server Native Client 10.0: Invalid object name 'gismo.RFCGIS.aocondo'.] [gismo.RFCGIS.aocondo][STATE_ID = 19653]
Failed to execute (TableToTable).

I use sdetable -o delete to remove the tables and that, too, has been working fine, so I don't think there's a problem with SDE. I've also checked to see if there are any remnants of the table names in sde_table_registry, sde_layers, sde_column_registry, and even gdb_items (though the tables should not have been registered with the GDB). I also checked against various owners (dbo, rfcgis, and sde) to see if it was a problem with the file existing but owned by a different user and that's not it either. Oh, and it does work when copying from DBF to DBF, so it's not contention on the source file. I've run out of ideas and hope someone can suggest another direction to take with this.


texts: How do I get contours from a scanned map?

texts: I have a scanned map which is a .tif file. I have geo-referenced it, and now I want to render it in 2.5D. Can i do this in QGIS? The ENC file is not available for the map I am using. The map already has contour lines on it, it is a colored sea map with black lines for contours. Land is denoted in yellow and green. Unfortunately due to company policy I cannot upload it.


texts: How to scale SVG symbol based on zoom level in QGIS?

texts: I am wondering if there is a way to scale SVG files I have stored as blob in a Spatialite database according to the current zoom level in QGIS. My hope is to have the SVG icon resize itself when the map is zoomed in/out. I understand there is a size option in the data defined option, but it does not seem to provide a dynamic sizing option. I am open to options.


texts: AttributeError: Object: Environment  not found from ArcGIS 9.3 Geoprocessor?

texts:
This problem is occurring in ArcGIS Server 9.3, Python 2.5.2

I am exposing a typical python script as a geoprocessing service.  In this script, I frequently use the gp.ListFeatureClasses() method, which searches the current workspace as shown in snippet below.

import arcgisscripting
gp = arcgisscripting.create(9.3)
gp.workspace = "C:/someworkspace/mygdb.gdb"
featureClasses = gp.ListFeatureClasses()

In the initial portion of the script, everything works as expected.  But in later portions of the script, the service crashes when I call gp.workspace reporting the following error:

  AttributeError: Object: Environment <workspace> not found

My initial guess was that maybe the workspace was being set to a non-existent workspace or gdb, but I have been able to verify that the workspace is always infact valid.

Digging around a bit has shown that others have experienced the same problem.  However, none have been able to set the workspace AT ALL.

Has anyone ever experienced such issues?


texts: Using UpdateCursor based on current selection of features

texts:
What I want to do is to manually select two features with the mouse (e.g. point or polygon), each in separate geodatabases to then copy the attributes values from one to the other (for only these two selected objects).

So far I have only managed to copy all attributes from one to the other.  I am stuck at how to make update cursor run only on selected objects

I am trying to use something like this to get the selected objects, but this selects the entire layer:

mxd = arcpy.mapping.MapDocument("CURRENT")
layers = arcpy.mapping.ListLayers(mxd)
sel_layers=[str(i.name) for i in layers if arcpy.Describe(i.name).fidSet]

And my update cursor is this:

fc_search = r'C:\PythonFolder\Othergeodatabas.gdb\Polyline'
fc_destination = r'C:\PythonFolder\LKpolygon.gdb\Avloppsbrunn'

with arcpy.da.SearchCursor(fc_search,["DIM","MATERIAL"]) as search_cur:
    for search_row in search_cur:
        with arcpy.da.UpdateCursor(fc_destination,["Tillkomst_plan","Notering"]) as upd_cur:
            for upd_row in upd_cur:
                    upd_row[0] = search_cur[0]
                    upd_row[1] = search_cur[1]
                    upd_cur.updateRow(upd_row)

This code works, but it applies changes to attributes for the entire layer, while I would like to only change attribute values for selected objects.


texts:
I'm new to QGIS and have no programming skills.

Is there an easy/uncomplicated way to create multiple buffers around a line feature?


texts: Find Closest Route from single source to several different destinations in PostgreSQL

texts:
I've made network nodes in PostgreSQL SQL based on Osm2po. I can find the closest way among two points based on the pgr_dikstra algorithm but I want to find the closest route between one single source point and several destination points just in one query and choose the closest one. How can I implement this query? For example, suppose I want to find the closest vehicle to me when there are many available vehicles

This is my code

SELECT * FROM pgr_dijkstra(
        'SELECT gid as id, source, target, length_m AS cost FROM ' || $1,
        -- source
        (SELECT id FROM ways_vertices_pgr
            ORDER BY the_geom <-> ST_SetSRID(ST_Point(x1,y1),4326) LIMIT 1),
        -- target
        (SELECT id FROM ways_vertices_pgr
            ORDER BY the_geom <-> ST_SetSRID(ST_Point(x2,y2),4326) LIMIT 1),
        false)

but this code works just for two points


texts: Fail to display correctly newly inserted data in PostGIS on Geoserver

texts:
I have created a postGis store in Geoserver. I have created layers with points and lines and I display them on an OpenLayers based application.
I also use an FME workflow to insert new data to my postGis table. I expected that my geoserver layers would be updated automatically but it didn't happen.

The problem was that I had to set new boundaries in the layers properties. After this I could see some (not all) of my newly inserted data on the map.
In different zoom levels I had different information displayed. For example in the upper zoom lever I could see the whole information but zooming in I was seeing less (see figures).

I thought it might be related with caching. I have used the geoserver interface to seed-generate missing tiles and according to geoserver the process was successful.

Lastly I have restarted geoserver and tomcat7 service without any success either.

When I select layer preview in


texts: Weighted kernel density estimate mapping in R

texts:
I am trying to reproduce ArcGIS' kernel density in R. Specifically, I'd like to make a weighted kernel density map, like you see in Census OnTheMap (example). I am running into a lot of dead ends, where you can easily make an unweighted density map, or you can easily plot a weighted density, but it's difficult to make a weighted density map.

The closest examples use smoothing over a specified area or expect you to already have a raster set up. I'm new to raster, so it's difficult to make sense.

I do know that ArcGIS defaults the bandwidth to Silverman's rule of thumb (bw.nrd0).

So how do we take points (say polygon centroids) with a value and make a map plotting the density of those values?


texts: filling the attributes according to intersection of two layers (Postgresql)

texts:
How to automatically fill attributes of a point layer according to the exact position in the polygon layer?

I mean how to write a command in PostgreSQL to automatically fill the attributes by the polygon layer (corine land cover) where is the point located?


texts: Calculating dominant angle of a polygon in QGIS without main_angle function

texts:
How can I calculate dominant angle of a polygon in QGIS without main_angle function?
I have a shapefile of Oriented Minimum Bounding Box rectangles (accessible from the following link: https://drive.google.com/drive/folders/1Xr_0Apg9s7s0daevIUctWTr0-WCABC83?usp=sharing) produced by using as Input Layer a shapefile of irregular polygons (accessible from the following link: https://drive.google.com/drive/folders/1BKp5VbR1aMAzCObGprvqeL0WSMqpTX32?usp=sharing)

Figure 1: Shapefile of Oriented Minimum Bounding Box

Figure 2: Shapefile of Irregular Polygons
I need to find the main_angle of every polygon in the Oriented Minimum Bounding Box layer.
I used the main_angle function available in QGIS 3.16 onwards.
main_angle function returns the main angle of a geometry.
However, when I constructed centerlines using the centroid of the polygons (in this case rectangles) and main_angle values computed using the main_angle
function in QGIS with the following expression:

intersection (
extend (
    make_line (
        centroid ($geometry),
        project (
            centroid ($geometry),
            25000,
            radians("main_angle")
        )
    ),
    25000,
    0
),
$geometry

)
Main_angle of the polygons computed with main_angle function is not at all times
the angle of longest collection of segments as in the image below encircled below:

What expression can be used such that main_angle is always the angle of longest segment of the Oriented Minimum Bounding Box in a project whereby main angle of every Oriented Minimum Bounding Box is the determining factor of accuracy of end results.
I found this comment from the following link: https://georezo.net/forum/viewtopic.php?id=119435
You need an expression dependent on the height and width data:

The end result I am seeking is the overall orientation of all ploygons collectively as illustrated in the image below from the Line Direction Tool QGIS plugin. By unfortunately provides no values.

I used XLSTAT to use Principal Component Analysis but did not achieve conculsive results. Computed results and report is accessible from following link: https://drive.google.com/drive/folders/1CHDfJ9poD8tk-KQLmRV8NbilZCyGtTs5?usp=sharing


texts: How to retrieve raster's metadata from PostGIS, using Python?

texts:
I am working on a PostGIS database, and I would like to store GeoTIFF raster images in some of my tables. Currently, I use GDAL to convert / read my images, and psycopg2 to connect to a SQL server. However, when I want to get a raster from the database, the metadata is missing.
I am using ST_FromGDALRaster to store raster data to the database (in bytes format), and ST_AsGDALRaster(my_tile, 'GTiff') to retrieve it.
Here is a code snippet that illustrates my issue:
First, let's create a naive gdal.Dataset object using numpy.
import numpy as np
from osgeo import gdal, osr

def create_ds():
    # GeoTIFF data (N_band, N_x, N_y)
    data = np.ones((3, 1000, 1000))  # naive data
    epsg = 3857  # mercator
    # Coordinates
    X = [-10, 10]
    Y = [-10, 10]
    # bbox
    Nb, Nx, Ny = data.shape
    xmin, ymin, xmax, ymax = np.min(X), np.min(Y), np.max(X), np.max(Y)
    xres = (xmax - xmin) / float(Nx)
    yres = (ymax - ymin) / float(Ny)
    geotransform = (xmin, xres, 0, ymax, 0, -yres)
    # Temp file, so nothing is saved on the disk
    vsimem = '/vsimem/temp'
    # Create a gdal.Dataset
    ds = gdal.GetDriverByName('GTiff').Create(vsimem, Ny, Nx, Nb, gdal.GDT_Float32)
    ds.SetGeoTransform(geotransform)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    ds.SetProjection(srs.ExportToWkt())
    # Write data for all bands
    for band in range(Nb):
        ds.GetRasterBand(band + 1).WriteArray(data[band])

    return ds

# Create the dataset

ds = create_ds()

# Modify the metadata
ds.SetMetadataItem('NewMetaKey', 'NewValue')
ds.SetMetadataItem('NewMetaKey2', 'NewValue2')

ds.GetMetadata()

# Out: {'NewMetaKey': 'NewValue', 'NewMetaKey2': 'NewValue2'}

Then, let's implement a function to convert gdal.Dataset to bytes, so it will be easier to store in PostGIS:
def to_bytes(ds):
    # use temp files
    vsimem = '/vsimem/temp'
    gdal.GetDriverByName("GTiff").CreateCopy(vsimem, ds)
    f = gdal.VSIFOpenL(vsimem, 'rb')
    gdal.VSIFSeekL(f, 0, 2)  # seek to end
    size = gdal.VSIFTellL(f)
    gdal.VSIFSeekL(f, 0, 0)  # seek to beginning
    dbytes = gdal.VSIFReadL(1, size, f)
    gdal.Unlink(vsimem)
    return dbytes

# Convert the dataset to bytes

dbytes = to_bytes(ds)

Then, convert bytes data to a gdal.Dataset:
def from_bytes(dbytes):
    vsimem = '/vsimem/temp'
    gdal.FileFromMemBuffer(vsimem, bytes(dbytes))
    ds = gdal.Open(vsimem)
    gdal.Unlink(vsimem)
    return ds

# Convert bytes to a dataset

ds2 = from_bytes(dbytes)

# Get the metadata
ds2.GetMetadata()

# Out: {'AREA_OR_POINT': 'Area', 'NewMetaKey': 'NewValue', 'NewMetaKey2': 'NewValue2'}

Next, add the data to a PostGIS table, named 'table':
import psycopg2 as psql

with psql.connect(**config) as conn: # To autocommit / rollback
    with conn.cursor() as curs:
        curs.execute("INSERT INTO table(my_tile) VALUES (ST_FromGDALRaster(%(byte)s))", {'byte': dbytes})

Then, retrieve the data:
conn = psql.connect(**config)
curs = conn.cursor()
curs.execute("SELECT ST_AsGDALRaster(my_tile, 'GTiff') FROM table;")

ds3 = from_bytes(bytes(curs.fetchone()[0]))

# Get the metadata
ds3.GetMetadata()

# Out: {'AREA_OR_POINT': 'Area'}

And at the end, the metadata does not match the initial values.
How can I modify my code so the metadata is not changed ?


texts: Topological Editing of Two different layers in QGIS

texts:
Is there a way to apply topological editing to move vertices between two coincident features that are not in the same layer?

For example, if I have a crop land polygon coincident with a soil polygon (from a crop layer and a soil layer) and I want to reshape the border between the two using the node tool, is that possible without moving a vertice in crop, changing to editing soil, and moving the vertice to be coincident with crop?

What if there is a line road coincident with the edge of a polygon (i.e. a polyline layer rather than a polygon layer), again can I reshape the two at once?

I am aware of the topological editing, but that appears to only work between two polygons of the same layer.


texts: Updating GeoJSON of countries with new and different property for each country

texts:
I am trying to add a new property to my countries' GeoJSON.

    var countries = {
    "type": "FeatureCollection",
    "name": "countries",
    "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
    "features": [
    { "type": "Feature", "properties": { "NAME": "Fiji", "POP_EST": 920938 }, "geometry": { "type": "MultiPolygon", "coordinates":
{ "type": "Feature", "properties": { "NAME": "Tanzania", "POP_EST": 53950935 }, "geometry": { "type": "MultiPolygon", "coordinates":
{ "type": "Feature", "properties": { "NAME": "W. Sahara", "POP_EST": 603253 }, "geometry": { "type": "MultiPolygon", "coordinates":

The property would be different for each country, EG 10 cases for Fiji, 20 for Tanzania.

I have tried

countries.properties.Fiji.cases = 10;


texts: Monte Carlo simulation using QGIS and pgRouting on optimal sidewalk construction

texts:
I am new to spatial analysis and would appreciate some general direction on a project I am attempting, outlined below (I am starting from scratch).

GOAL: To find the best locations to install 2000 feet of sidewalks in my hometown in order to connect the most households to the Central Business District (CBD), where "connected" means within 1.2 miles walking of the CBD. I have shapefiles showing existing structures (households), roads, and sidewalks (already installed).

Here is my proposed solution / thought-process:

Convert the in-place sidewalk network into a database of nodes
that are connected by weights (i.e. distances).  Is there a way to
directly do this in QGIS (or other program) by clicking on all
intersections?
Calculate the number of households that are within 1.2 miles
walking of the Central Business District (e.g. a lat-long point or
polygon) using the routing capabilities of pgRouting or something
else.  This will be the base case "household access" value.
Using the road layer as a guide, randomly place an additional
2000 feet (say, in 10 foot segments) of sidewalks onto the sidewalk
layer.  This is the equivalent of constructing a bunch of new
sidewalks arbitrarily.
Re-calculate the nodes and weights using the new pedestrian
network as in (1), and then re-calculate the number of households
that are now within 1.2 miles of the CBD as in (2).  It should
increase with the additional sidewalks.  Save the locations of the
additional sidewalks and the associated "household access" value to a file
(e.g. spreadsheet).
Repeat steps (3) and (4) 10000 times, similar to a Monte Carlo
simulation.  Using the 10000 sets of data points, choose the sidewalk placement locations that maximize the number of households within 1.2 miles of the CBD.

Does this thought process sound realistic? Does anyone have any suggestions?

-- I would like to accomplish this using some combination of QGIS and R, however I am open to learning PostGIS and/or Python (or anything else) to achieve the goal.


texts:
Based on a previous post, I tried to implement this in the expressions tab:

CASE
  WHEN "dBM" <= 0 THEN '255,0,0'
  WHEN "dBM" > 0 AND "dBM" <= 20 THEN '255,255,0'
  ELSE '0,255,0'

However, I get an error like this:

Parser Error: syntax error, unexpected $end Eval Error: No root node! Parsing failed?

Desperately need this to work...


texts: Extracting raster cells from imagery based on RGB value using ArcGIS Desktop?

texts:
I've got some aerial photography (format: IMAGINE Image) which is detailed enough to see to fire hydrants and a shapefile of fire hydrants which was digitized based on some CAD data.  My project is to QC the fire hydrant shapefile locations by spot checking with the imagery. Is it possible to extract RGB values from the imagery that match the fire hydrants RGB values?

Fire hydrant (red):

I'm looking to select cell values (which i have written down) that are similar to the cell values of the fire hydrant and then extract those cells (either a raster file of those cells or a shapefile which plots location of similar cells). I'm assuming that there are very few features in the imagery that have similar reddish RGB Values.

After looking at Arron's Answer and playing around with both the supervised and unsupervised approach (see my comment below), i was unable to get the tool to perform exactly what i wanted until I started too look at the confidence raster output from the Maximum Likelihood Classification tool.  I am not to sure what exactly the raster output is supposed to symbolize but from just looking at the level 14 cell value, but it captured the all the fire hydrants.

raster output:

location of fire hydrants:

My next step is to use the tool Raster to Polygon convert the fire hydrant raster footprint to a polygon.  I am running the tool right now and in the past hour it has only made it to 11%.  The area I am working in is large, at 1x1 mile so I understand that it could take some time if there is a lot of small raster cells to convert to a polygon.  Are there any suggestions on running some tools to clean up some of the data so it speeds up the process of converting the raster to a polygon? I may not have a need to run this tool this since the raster footprint did such a good job at capturing the fire hydrants but I am interested in speeding up the process since I foresee this being used in other applications for the future.


texts: Exporting lines of .lyr file using ArcPy?

texts:
After I finished make an OD cost matrix layer with function

outNALayer = arcpy.na.MakeODCostMatrixLayer(inNetworkDataset, outNALayerName, impedanceAttribute, 15, "", accumulateAttributeName)

the result was saved as an .lyr file. But I want to get the lines result which is corresponding to the object circled in the picture.

(in Arcmap)

Is there any function I can use or are there any other method to save the OD matrix result?


texts: How to extend a custom NASolver in ArcGIS Server 10.1?

texts:
The past year we developed a custom NASolver with C# , it works fine and that job has finished.

Now we are planning to use it from ArcGIS Server. The problem is that we have no information for that. The creation of a new Server Object Extensions could be a solution, but it will take time...

We want to know if exists another way to do that by extending the iNAServer object.

Thanks in advance for your help.


texts: Plot the nodes of the shortest path found by dijkstra algorithm

texts: I am using qgis and pgrouting plugin to find the shortest path between two nodes and display the route. Currently, the plot shows the shortest path. However, I also want to know which nodes the path passes through. So is it possible to actually get this information?


texts: Interpolation of point data with different scope

texts: I need to interpolate a raster from point data (air pollution stations data). What I want to achieve is to distinguish the station's scope (or range or scale). For example some station's data are valid for area with radius of aprox. 1 km around the station and for some other one it is 10 km radius. I have this in separate attribute. I am looking for a way how to make this affect the interpolation results. I am able to work in ArcGIS or GRASS on this task. I am familiar with IDW, Kriging and some other interpolation methods but I couldn't solve this so far.


texts: Preventing symbology reclassification from adding 0s to Range in ArcGIS Desktop?

texts:
How do I prevent symbology's "Natural Breaks" from automatically changing my range to include multiple 0s after my decimals?

Example: Range was originally 1 - 249. I switched to Natural Breaks and the range is now 1.0000000000 - 250.00000, then 250.0000001 - x. It should be 251.

In Range my counts now have just enough zeros to make both columns have the exact same number of digits and line up perfectly.


texts:
I have a map of a region and I am plotting pie-charts (as diagrams) over every country in that region in QGIS (1.9.0). I also need to connect those pie-charts with lines (which I am creating in the Mapping Toolbox of MATLAB as a separate shapefile).

Can you please tell me if there is a way of knowing in advance what coordinates of a country QGIS shows its pie-chart at? Also, can I capture those coordinates automatically (even before the pie-charts are drawn) just as I would normally do using the Coordinate Capture plugin? That way, I can create the lines from one point to another before even plotting the pie-charts.


texts: Interpolating shapefile of points to grid using Java?

texts:
I am looking for a method/prog/library(java) to interpolate points from a shapefile to a grid (file). I have seen arcgis can do that but i am looking for something free.

More specific: i have shapefiles(256 or more files) with points, those points have a value of temperature. Now i want to make a file with the values of a grid converted out of this points. The interpolation method is not very important for this project. The grid would be something of 10x10.

I would like to do this in java and iterate over the files and convert them to grids. I presume there is a java library that can do this.


texts: Error in raster histogram in QGIS3.0 with PyQGIS

texts:
I am working on a plugin, one of its functions is to create the raster layer histogram.
I generate the histogram using the histogram method of QgsRasterInterface and graph it well, without problems in QGis3.4.
But in QGis3.0 it does not allow to generate the histogram.
This is a recognized error? when creating the histogram in QGis3.0?
capa=iface.activeLayer()
provider=capa.dataProvider()
provider.initHistogram(QgsRasterHistogram(), 1, 20)
histogram = provider.histogram(1,20)
y=histogram.histogramVector

Details:

If I get correctly the statistics with bandStatistics
The values of histogram.minimum and histogram.maximum do not coincide with the maximum and minimum of the statistics.
Try assigning the maximum, minimum and extension parameters. As well as placing None,None, but I still did not get a result
I try to create the histogram without initializing it or defining no data, could this be the error in QGIS3.0?


texts: How can I do spatial queries in client side?

texts:
My current project is Map visualization. In this application I am taking the Input as JSONArray/CSV file and the outputs are Heat Map, pin Map, Thematic Map etc.. So, here I am using postGIS for only to store boundaries(like state/county/city..). Also the data from the user will remain in his place, that means I am not storing the user input data in database.

How can I do some spatial queries/functions in client side?
ie, functions like st_buffer, st_intersect, st_contains, st_distance, st_equals, st_    boundary etc....

I am using Openlayers 2.11, Geoserver, postGIS


texts: Creating definition query on layer excluding decimal values in ArcGIS for Desktop?

texts:
I am working on a layer of spot height in ArcGIS 10.3. The attribute table of this layer contains a field called 'Elevation' wich one provide the altitude of these spots height. These values can be decimal (with a comma) or integer.

What I want is make a definition query in view to show spots height with integer altitude only. I don't want to represent spots height with a comma in their elevation.

I have tried seleral query expression in the layer properties of my layer but it doesn't run like this one : ALTITUDE <>%.%


texts: Determining the camera direction of a geotagged image

texts:
Is there an easy way to determine the camera direction of a geotagged image?
Maybe there are several methods for this.


texts: Arcgis 10.2.2 soap API Error {401 - Unauthorized}

texts:
I'm trying to consume secured map service via arcgis 10.2.2 soap api. I've configured my web adapter to negotiate windows domain security and everything looks fine when I check over the rest api. However I get a this exception when I try to consume the same service with soap api. I'm getting 401-unauthorized at getDefaultMapName method.

 mapServer = new MapServerBindingStub(uri, username, password);
 defaultMapName = mapServer.getDefaultMapName(); //<error here 401-unauthorized>

Error

Caused by: com.esri.arcgisws.runtime.exception.ArcGISWebServiceException: HTTP s
tatus code {401 - Unauthorized} received from server
        at com.esri.arcgisws.runtime.transport.http.HttpClientConnectionChannelI


texts:
Does anyone know where I can find the converter from .map to .sld?

I have a big mapfile of styles for OSM and I have the styles mapnik.xml for OSM but I use Geoserver for my work. Or does anyone have some sld styles of OSM for PostGIS layers planet_osm_point and planet_osm_polygon? (I already have the style for planet_osm_line.)


texts: How does ArcGIS Handle Multiple Instances of the Same Attachment Across Different Features?

texts:
I am considering enabling attachments in several different feature classes.  Our Enterprise Geodatabase represents as-built pipeline system features, such as joints of pipe, welds, and valves.

What I would like to do is store pictures on some features and PDFs (such as MTRs) on others.  However, I am concerned that the nature of these features will contribute to unnecessary bloat in the size of the database.  I will frequently need to store the same single file against multiple different features.

For instance, if I have a valve set of 5 valves all pictured in the same single photo, will using the attachments capabilities result in this single photo being stored 5 times?

Or if I have a batch of 75 joints of pipe, all with the same heat number, and thus all with the same MTR pdf file, will that pdf file be stored 75 times?