texts: Overpass API: Loop through cities and find total length of road centerline per city

texts:
I am trying to prepare a Overpass QL such a way that it iterate through a list of City Names and find sum of lengths of all the road centerline features which are falling inside the city administrative boundary.
The idea is not to download to entire data rather to get only required total road length per city wise.

New York City, New York, USA
Toronto, Ontario, Canada
Jersey city, New Jersey, USA
<City_Name>, <STATE/PROVINCE>,

Above is the example input and below is sample QL with BBOX but need to automate either with overpy Python AOI or writing single overpass QL for entire list of cities.
[out:json][timeout:25][bbox:{{bbox}}];
(way["highway"="primary"];
 way["highway"="secondary"];
 way["highway"="tertiary"];
 way["highway"="unclassified"];
 way["highway"="residential"];);
(._;>;);
out;

Can someone give me pointers?


texts: How to create a file with the names of cities and the state of the city from OpenStreetMap?

texts:
Is it posible to extract data from OpenStreetMap that only contains the name of a city and the name of the state from these city?

I want to build a CSV file, where each row contains "city;state".
For example:

San Francisco;California
Annapolis;Marryland
....

At the moment I'm using osmfilter and osmconvert to create a file from all cities.

osmfilter file.o5m --keep="place=city or ( place=town and population>=1000 )" --drop-relations --drop-ways > cities.osm
osmconvert cities.osm --all-to-nodes --csv="name" --csv-headline > cities.csv

But how can I add the information about the state / country?


texts: How to enable a keyboard shortcut for a plugin

texts:
I've created a QGIS plugin using plugin builder which works just fine when I enable the plugin dialog from the toolbar button or the plugin menu item. The plugin has one dialog with a few combo boxes, and an ok/cancel button. The user makes combobox layer selections, then makes feature selections on the canvas, then presses ok to run some code based on those choices and selections.

However I can't seem to get a keyboard short cut enabled, one that would allow the user to keep the plugin dialog open and a keyboard shortcut that would accomplish the same thing as clicking the OK button but without closing the dialog.
I've had a look at this method here, but I couldn't get it to trigger the run function where most of my code sits. My initGUI function looks like this:

self.keyAction = QAction(u"Copy the Attributes", self.iface.mainWindow())
self.iface.registerMainWindowAction(self.keyAction, "F6")
self.keyAction.triggered.connect(self.run)

But that doesn't  seem to do anything when I press F6. I've tried setting the shortcut in Settings -> Configure Shortcuts, but that doesn't result in anything happening either. My unload function is also probably wrong as well:

def unload(self):
    """Removes the plugin menu item and icon from QGIS GUI."""
    for action in self.actions:
        self.iface.removePluginMenu(
            self.tr(u'&Copy Attributes'),
            action)
        self.iface.removeToolBarIcon(action)
    # remove the toolbar
    del self.toolbar
    #self.iface.unregisterMainWindowAction(self.keyAction)

I'm struggling to understand if I'm even going in the right direction with this registering and unregistering, if I want a single keyboard shortcut key to mimic the press of the OK button in my plugin dialog - while keeping the dialog box open. Keeping the dialog open and using a shortcut key would allow the user to perform multiple plugin actions on the same layers, but with different feature selections on the canvas.

Any thoughts on where my approach goes wrong?

Update: As suggested below, I've added a push button in the dialog which allows me to just run the code needed without closing the dialog:

self.dlg.pushButton.clicked.connect(self.run)

But is it still correct to try to use a shortcut using

self.keyAction = QAction(u"Copy the Attributes", self.iface.mainWindow())
self.keyAction.triggered.connect(self.run)
self.iface.registerMainWindowAction(self.keyAction, "F6")

...when I'm not using the keyboard shortcut to open the dialog, but rather run the code under self.run ? Because it doesn't work when I try this. I'm missing something obvious here.


texts: Calculate field in table A using multiple values from table B in ArcMap?

texts:
I currently have a point shapefile, and a table.  The tables are set up in a many-to-one relation with "PlotName" as the primary and foreign keys.

The shapefile contains the primary key (i.e. A West H1), the table contains many fields with the same key.

In the table, each row has both a key (PlotName), a value (PercentCover), and whether this value is Native or Invasive (0 or 1).

I would like to create a field in the shapefile (TotalPercentNative) that contains the sum of (PercentCover) based on whether it is native or invasive.

Do you have any suggestions?


texts: Calculating distance between linestring and polygon with geopandas


texts: Replace the attribute value of all shapes with specific value in QGIS

texts:
I'm pretty green when it comes to QGIS. I have a shapefile with a bunch of fields and farm codes of whats growing in them. However there are like 600 of them, so I want to go through and simplify. So  All values that represent a type of fruit tree will just be one value.
How can I implement a function along the lines of 'IF 'field value' = x,y,z THEN = Q?'
and where would I implement that? Somewhere in the attribute table right?


texts: Calculate elevation profile for a line from dem data

texts:
We are trying to implement a service which takes a line geometry as its input and calculates elevation profile along that line from dem data (img format).

How could we achieve this using postgis or arcgis stack?


texts: ArcMap 10.1 Spatial Join: 'Contains' and 'Within' question

texts:
I have two polygon datasets: electoral districts and postal code areas. I'm performing a spatial join to see which postal code areas are completely within a district (the area can touch the borders, but doesn't cross them; i.e. not situated in 2+ districts).

The way I see it I can do it two ways:

(1) Districts as target features, postal code areas as join features with a Contains Clementini match option

(2) Postal code areas as target features, districts as join features with a Within Clementini match option

I assume these are the same thing. Yet I have ran both and they yield different results.

Theoretically: why am I getting different results? And which way is deemed "more correct"?
Thanks for your input in advance!


texts: QGIS 3.4.3 Color region bounded by certain lines

texts:
How do I color in the certain regions that are bounded by lines?

I'm completely new to QGIS.


texts: Does OSM road segment has GPS information about the Intersecting roads in them?

texts:
I have been reading about OSM Road from accepted answers here like,

how-to-split-osm-roads-into-individual-segments-at-intersections

finding-all-segments-that-makeup-a-road-in-osm

and my question,

does-osm-road-segments-contain-intersections

The answers fell short in explaining whether the segments in an OSM Road belong to a single road or grouped together from multiple roads. Since OSM Roads have intersections(answered in third link above), does it mean OSM Roads actually contain roads where travelling in multiple directions is possible.

To put it clearly,

Is the road uni directional - one way or two way
Does the OSM roads has structures like T-Joints (where the road has information(GPS points) on roads connected to an intersection.). Usually it is not possible for us to travel in both directions at the same time. So having OSM ID for GPS points in different road segments where you cannot travel, does not make sense to me.
Where can I find the reference or documentation for these latent  information for OSM.

Note: If the question seems too abstract, let me explain, I am trying to understand the OSM, since I'm using it for my research.


texts: How to aggregate a list of raster in Composite Bands arcgis

texts: I created a list of raster using ListRaster. However, composite band tool does not support list input. I want aggregate the list of raster into one raster, is there any other tool to do it in ArcGIS(arcpy)?


texts: Automatic shapefile input for pre-constructed QGIS workflows?

texts:
I currently have two models running:

Take shapefiles and break them down into separate shapefiles based on their class, and save them in a specific folder;
Select the broken-down shapefiles with the classes I want and merge them. Then, after they are merged, dissolve and save as a new shapefile.

These models work fine but the problem is that I have a lot of files to run these processes on. Each run takes roughly an hour and I'd like to just leave the computer running for a day so I dont have to keep checking back on the progress.

Is there a way such that QGIS can automatically input those files into the model and iterate through, other than just adding more arguments and parameters to my existing model?


texts: What is the function of GeoServer's ImageMosaic "imageindex" field?

texts:
I'm trying to use GeoServer's ImageMosaic plugin to serve up some folders of netCDF files with a particular naming convention. Rather than try to rig up some collection of regex-based PropertyCollectors, I'd like to manually populate the granule index.

I've noticed that ImageMosaic requires the index table to have an imageindex integer field, whose values seem to form a globally unique sequence (i.e., unique across all variables in the mosaic.)

What is the purpose of this field, and is there any requirement on how I populate it?


texts: Merging polylines together keeping the attributes of the largest one

texts:
I have a set of polygons with polylines running through them. I am trying to assign the polylines the polygon ID that they pass through the most. I have already split the polylines to the polygon boundaries, merged the individual shapefiles together keeping seperate features and joined them spatially to the polygons to give them the polygon IDs.

The issue I now face is to dissolve the features based on their polyline IDs while making the features keep the attributes of the section that is longest in a polygon. For example, if 60m of a polyline falls in polygon 1 and 40m of it falls in polygon 2, the rejoined line keeps the attributes of the section in polygon 1. There may be any number of features to merge.

Notable Condition: Some polylines may have more length outside of a polygon than inside. For those lines, I need the line to keep the attributes of the one that passes through a polygon. For example, if 5m of the line passes through polygon 1 and 95m of it pass through empty space, the merged attributes should reflect the line passing through the polygon

Here is the basic logic that I believe will work. I am just not sure how to code it into python or if there are tools that can accomplish what I am trying to do.

Select two line features. When one line has lineid = 1 and the second has lineid = 1, check for polygon ids in both, if only one has a polygon id, merge based on the one with an id, if they both have polygon ids, compare the lengths and merge based on the one that is longest. Move on to the next feature.

I am not sure if they need to be compared one at a time or if the dissolve tool or something else incorporates this kind of code. Thoughts?

I am using ArcGIS 10.0


texts: QGIS Mosaicing with null values in images

texts:
I have a set of landuse landcover data which needs to be mosaiced.

Whenever i try using QGIS Raster>miscellenous>merge tool. giving no datavalue as '0', it returns either

(A)

 ' Traceback (most recent call last):
File "C:\PROGRA~2\QGISVA~1\bin\gdal_merge.py", line 509, in
sys.exit(main())
File "C:\PROGRA~2\QGISVA~1\bin\gdal_merge.py", line 495, in main
fi.copy_into( t_fh, band, band, nodata )
File "C:\PROGRA~2\QGISVA~1\bin\gdal_merge.py", line 241, in copy_into
nodata_arg )
File "C:\PROGRA~2\QGISVA~1\bin\gdal_merge.py", line 61, in raster_copy
nodata )
File "C:\PROGRA~2\QGISVA~1\bin\gdal_merge.py", line 101, in raster_copy_with_nodata
to_write = Numeric.choose( nodata_test, (data_src, data_dst) )
File "C:\PROGRA~2\QGISVA~1\apps\Python27\lib\site-packages\numpy\core\fromnumeric.py", line 296, in choose
return choose(choices, out=out, mode=mode)
MemoryError

or

(B)

' The process crashed sometime after starting successfully'

I am using QGIS 2.2 on Windows 8 PC.

I have a set of 30 images in .tiff format which needs to be mosaiced. I also tried doing it with 2-3 images but the same problem persists.


texts:
EDIT 10/20/2011:

The problem I was having below was caused by trying to install Python(x,y) after ArcGIS 10. I got a new machine yesterday and instructed IT to install Python(x,y) 2.6.5.3 (has correct numpy version for Arc). I had them run a non-scripted installation of ArcGIS 10 so we could make sure the installer didn't try to give us a second version of Python. I now have a single Python directory and ArcPy works (so far); as of right now I'm considering this issue resolved.

I have what would seem to be an issue with a clear and easy answer, but so far I haven't been able to find a suitable one.

Here is the situation:

I need to have both ArcGIS 10 and Python(x,y) installed on a server that is accessed by multiple users in order to accommodate everybody's needs. Because both of these pieces of software install their own version of Python 2.6, the libraries installed by one cannot be accessed by the other. In other words, if we install ArcGIS before Python(x,y) we can use all of the goodies like MatPlotLib, numPy, etc., but the arcPy module is unavailable and vice versa.

To solve this problem I tried adding the path to the various libraries to the python path of the "active" version of Python to no effect, ala this thread: Problems setting up eclipse (with PyDev) for ArcGIS 10.

We could possibly resort to installing ArcGIS 10 and then try to manually install all of the libraries that are included with Python(x,y), but this would not be an ideal fix. Alternatively we could install ArcGIS 9.3.x since that comes with Python 2.5, but I would prefer to have all of the functionality that is available with the new arcPy module. If someone out there has already encountered this and has any other suggestions I'm all ears!

We are working in the Windows Server 2008 R2 environment.


texts:
I'm attempting to make a flood plain map for a school project, and am using FEMA Flood Hazard Symbology for the legend. I've dissolved by a Flood Zone Field, but am left with two attributes both named "X" they are slightly different, but for what I need it does not matter.

I want to consolidate these attributes into one. So far I've tried to merge them using the editor toolbar. The process starts to run, but ArcMap locks up and eventually says the merge cannot be completed: Not enough memory. I'm running Arc on a computer with 8Gb or RAM. Is there another way to consolidate these into one, or a way to not make it lock up?


texts: how to combine chloropleth, bubble and cluster designs in CARTO?

texts:
To be concrete, I have a map of election results by precinct. As I do not have shapefiles for the precincts, I geocoded the polling places, but many precincts share polling places. I have the following design in mind, but as it goes beyond the available wizards in Carto, some guidelines about what is feasible with CartoCSS (and mabye how) would be welcome.

It would be nice if the precincts did not overlap. Some design where they are clustered by polling place or "jiggled" to appear around the polling place but not on top of each other would be nice, if possible.
The size of the markers should vary be the size of the precinct, just like in a bubble design.
The color/filling of the markers should vary by the vote share in the precinct, basically as in a chloropleth design.

+1. It would be nice if the map could be interactive and the user could toggle between two different vote shares coloring the same precincts with the same bubble sizes. (Basically, "No" votes in a referendum vs. invalid ballots.)

+2. A Torque-like design could replace point 3 (but still combined with points 1 and 2), where precincts bubble up over time according to their results. I was surprised this was not easier: the election results being numbers between 0 and 0.5, they did not automatically span the time of the animation, I am not sure if this is possible without manipulating the underlying data into a different range of numbers.

I hope at least points 1-3. are possible to combine.


texts: ArcMap10: Changing transparency of a polygon affects symbol in table of contents but not in the map

texts:
I am working with polygon shapefiles both that I have created and that colleagues have sent me. When I try to change the transparency of the polygon through Layer Properties > Display >"Transparent: __%", the symbol in the Table of Contents on the left is changed to be more transparent, but the shapefile in the map view still has 0% transparency.

I also noticed in Layer Properties > Symbology > "Advanced: Transparency", you can "Vary feature transparency based on field values in percent." The shapefile that I created does not have much attribute data (just FID, Shape, ID), and the shapefiles from my colleagues have a few more attribute fields such as Shape_Length and Shape_Area.

However, when I toggle Symbology > "Advanced: Transparency" to different attribute fields, sometimes the transparency of the shapefile changes even though I did not change Display > "Transparent: __%". The transparency changes to some random (random to me, because I don't know why) amount. This random change is not consistent across attribute fields in different polygons.

I have tried changing the transparency through the Effects Toolbar, but I have the same problem where the symbol in the table of contents has an updated transparency but the polygon in the map still has 0% transparency.

Is there somewhere else transparency information is being stored/edited?

Your comments are appreciated! Thank you so much. I learned to use ArcMap 9 and am getting back into GIS through ArcMap 10.


texts: Dividing polygon layer into multiple layers by values in attribute table column using ArcGIS for Desktop?

texts:
I have polygon layer with multiple types of area (forest, grass, and many more), and I´d like to divide this layer into many layers by the area type.

I´ve tried Split function but I receive an error 000354 - because the type of areas are in Czech and contain national characters.

Please does anybody know how to do this??

For example how to automatically rename the types of areas to valid form, or how to write an python script which creates this layer for each of the area type in with definition query  for layer?

I'm using ArcGIS 10.1 for Desktop.


texts: Batch convert dgn files to shp files using ogr2ogr

texts:
We have hundreds of dgn files which have to be converted to shapefiles. I am trying to use these commands to batch convert all the dgn files in a folder.

for %%f in (*.dgn) do ogr2ogr -f "ESRI Shapefile" ".\point\%~nf.shp" "%%f" -where "OGR_GEOMETRY = 'POINT'" -skipfailures
for %%f in (*.dgn) do ogr2ogr -f "ESRI Shapefile" ".\linestring\%~nf.shp" "%%f" -where "OGR_GEOMETRY = 'LINESTRING'" -skipfailures
for %%f in (*.dgn) do ogr2ogr -f "ESRI Shapefile" ".\polygon\%~nf.shp" "%%f" -where "OGR_GEOMETRY = 'POLYGON'" -skipfailures

But I am getting an error of

  "|f" is not recognized as an internal or external command, operable
  program or batch file

Could anyone tell me what the problem with this command is?

Update: The problem is with the bat file format. It has to be saved in ANSI & not in other format's. The code given by Zoltan in the answer works perfectly

https://stackoverflow.com/questions/27205097/ansi-unicode-conflict-in-batch-file


texts: Google Earth Engine, Image Collection Statistics?

texts:
I've wrot below code for time series processing via modis LST product, but it shows this error:

ImageCollection (Error)
Error in map(ID=2018_06_02):
Image.constant: Parameter 'value' is required.

Google Earth Engine Code: https://code.earthengine.google.com/69da688a94c65e05e3cac421d5bbd938

Map.centerObject(table);
Map.addLayer(table);

// test for tci
var modis = ee.ImageCollection('MODIS/006/MOD11A2')
.filterBounds(table)
.filterDate('2018-01-01','2019-01-01')
.select("LST_Day_1km")
.map(function(img){
  return img.multiply(0.02)
  .copyProperties(img,['system:time_start','system:time_end']);
});

var tci = modis.map(function(img){

  var min = img.reduceRegion(ee.Reducer.min(),table,1000).get('LST_Day_1km');
  var max = img.reduceRegion(ee.Reducer.max(),table,1000).get('LST_Day_1km');

  var index = img.expression(
    '(max-LST)/(max-min)'
    ,{
      'LST': img,
      'min': ee.Number(min),
      'max': ee.Number(max),
      }).rename('TCI')
      .copyProperties(img,['system:time_start','system:time_end']);
  return index;
});

print(tci)


texts:
I have two shapefiles containing each one set of points, set A and set B. One set is my random selection (set A), the other is a list of communities (set B).

For each point of the random selection (set A) I want to find the nearest point pertaining to


texts: Getting the layer styles out of a *.QGS without having the layers itselves

texts:
Is it possible, to restore the layer-styles from a *.QGS project file?

I have unhandled layers in my project. As I'm starting the project, QGIS is asking me for the source of the layers but I don't have them anymore. So I ignore it and the project will be - of course - loaded without any layers and without their styles. But I just need the styles - how do I get to them as style file?


texts: Using field in place of value in WHERE clause of arcpy.da.SearchCursor()?

texts:
I am developing an ArcPy script in ArcGIS Desktop 10.7 and using a file geodatabase feature class. I have a field called facility_id. I have tried using this field as a string field or as an integer field but have been miserably unsuccessful at using this field as a value in the arcpy.da.SearchCursor function. For my latest attempt I generate the value in facility_id from an integer field called x.

x = 123456
facility_id = "'" + str(x) + "'"
print facility_id gives '123456'

I would like to do this:

cursor = arcpy.da.SearchCursor(facility_fc, ("SHAPE@XY"), "ASUFAC = '$facility_id'")

(ASUFAC is a string field)
When I execute this line of code, nothing is selected. If I execute the following line of code, one record is selected from the feature class, and is correct:

cursor = arcpy.da.SearchCursor(facility_fc, ("SHAPE@XY"), "ASUFAC = '123456'")

What is the correct syntax to use in order to have a value from the facility_id field be evaluated correctly in the where clause of the arcpy.da.SearchCursor function?

I've reviewed the thread "Using numeric variable in where clause of Arcpy SearchCursor" and adapted what was given there but it only generates an error: Invalid SQL statement was used


texts: Requirements for developing with ArcGIS for Android?

texts:
I want to develop a GIS application based on the ArcGIS for Android SDK.

What are the hardware and software needed by a client to run this project after building it?


texts: Clipping raster exactly to outline of shapefile using QGIS?

texts:
I have clipped the shape I want but it looks messy with ragged edges. I want to remove these and have clean edges.

I have tried doing this on QGIS but have come to the conclusion that it is not possible. Hopefully I am wrong. I have access to QGIS, GIMP photo editor, paint , R and Python. Unfortunately I don't have ArcGIS where I believe it is easily possible.


texts: Showing web page rather than just URL in ArcGIS Online pop-up box?

texts:
How can I show an actual web page (rather than just a hyperlink) in the pop-up box of ArcGIS Online?

In the ArcGIS Online "Pop-up Contents" section adding a url in the "Custom Attribute Display" just shows a link, and clicking on it opens the url in different page.

Is there any way to show the web page directly in the map pop-up like happens for images or charts when we click on a layer, as in the mockup below?


texts: Describing element in geodatabase using ArcPy

texts:
With help from people on this site, I have the following code to look for all of the element types in the geodatabase and put them into a report.

walk = arcpy.da.Walk(workspace, topdown=True)
    for dirpath, dirnames, filenames in walk:
        for filename in filenames:
            desc = arcpy.da.Describe(os.path.join(dirpath, filename))
            elements[desc['dataType']].append(desc)

report = ""

for element_type, element_list in elements.items():
    report += "{}:\n".format(element_type)
    for element_description in element_list:
        report +=
ntpath.basename("\t\t{}\n".format(element_description['featureClassType']))

Now I want to describe each of the element in detail. For example, if I want to include in the report what shape type (point/line/polygon/annotation) each feature class in the geodatabase is, I think I can use featureClassType, but I still do not understand how I can put it in to the code above.


texts: Will rebuilding a spatial index recalculate the feature extent?

texts:
ArcGIS 10.3.1

ArcSDE 10.2.2

SQL Server 2012

I'm reading documentation on how to Rebuild a spatial index. Will rebuilding the spatial index automatically recalculate the feature extent of the feature class?


texts:
It seems that in API for QGIS 3 VectorWriter was replaced:

from processing.tools.vector import VectorWriter returns
ImportError: cannot import name 'VectorWriter'.
How do I rewrite following things for PyQGIS 3?

writer = VectorWriter(output, None, fields,
                              QGis.WKBMultiPolygon, input_layer.crs())


texts: Autofilling polygon with semi-regular spaced points

texts:
I'm looking to fill a polygon with points with a regular spacing, but the points don't need to be in line vertically.
So, in the image below, the points are all spaced out evenly on an x/y spacing, but while the rows are all in line, they have been shifted left or right to match the polygon shape and fit in the maximum points on each horizontal row.

This was done on another piece of software, which we no longer have access to, but I'm looking to do the same in QGIS.
Any ideas?
So thanks to @Babel's help, I'm able to get somewhere, but where two polygons meet there is an overlap, the intersection expression worked for displaying the lines but not in the 'geometry by expression' output.
Id like to keep the lines within their own polygon, rather than the overlap.


texts:
Can anyone recommend a good reference or tutorial that explains the principles of spatial references, projections, and coordinate systems?

We are a small team of programmers working with GIS technology, and we often get snagged when working with data in different SRS/CS because we don't understand the concepts.


texts: Using simple features (sf) in R, how do I erase polygons overlapping with another layer

texts:
With a regular spatialDataFrame I could use erase(x, y) or x-y. How can I perform the equivalent using sf objects?

As an example:

##Create layers based on this answer: http://stackoverflow.com/questions/26620373/spatialpolygons-creating-a-set-of-polygons-in-r-from-coordinates
set.seed(145)
square1 <- t(replicate(50, {
  o <- runif(2)
  c(o, o + c(0, 0.1), o + 0.1, o + c(0.1, 0), o)
}))
ID <- paste0('sq', seq_len(nrow(square1)))

polys1 <- SpatialPolygons(mapply(function(poly, id) {
  xy <- matrix(poly, ncol=2, byrow=TRUE)
  Polygons(list(Polygon(xy)), ID=id)
}, split(square1, row(square1)), ID))

polys.df <- SpatialPolygonsDataFrame(polys1, data.frame(id=ID, row.names=ID))

square2 <- t(replicate(5, {
  o <- runif(2)
  c(o, o + c(0, 0.1), o + 0.1, o + c(0.1, 0), o)
}))
ID <- paste0('sq', seq_len(nrow(square2)))

polys2 <- SpatialPolygons(mapply(function(poly, id) {
  xy <- matrix(poly, ncol=2, byrow=TRUE)
  Polygons(list(Polygon(xy)), ID=id)
}, split(square2, row(square2)), ID))
polys.df2 <- SpatialPolygonsDataFrame(polys2, data.frame(id=ID, row.names=ID))

#plot the layers
plot(polys.df, col=rainbow(50, alpha=0.5))
plot(polys.df2, col="black", add = T)

#Compare methods - erase and difference
erasedPoly <- erase(polys.df, polys.df2)
plot(erasedPoly, col=rainbow(50, alpha=0.5))

polys.df_sf <- st_as_sf(polys.df)
polys.df_sf2 <- st_as_sf(polys.df2)
diffPoly <- st_difference(polys.df_sf, polys.df_sf2)
plot(diffPoly, col=rainbow(50, alpha=0.5), max.plot = 1)

Erase here does what I want, but st_difference doesn't. My workaround, if there is no in-built function in sf, will be to convert back to a Spatial PolygonDataFrame.


texts: Create 6 band IMG file with GDAL out of two JPG files, one with WorldFile coordinates

texts:
I have two JPG files

1 ortomap: http://imageshack.us/photo/my-images/703/ortofoto1.jpg/

with WorldFile:

0.683772574816354
0.000000000000000
0.000000000000000
-0.683739485024190
395057.114562241190000
741149.458682749890000

2 aerial photo in form of JPG: http://imageshack.us/photo/my-images/688/aerialphoto.jpg/

How can I create a 6 band IMG out of those two JPGs but according to the coordinates from the first JPG WorldFile using GDAL?

I am after some exemplary commadline commads


texts: Importing XML document to file geodatabase

texts:
I have never used xml data imports.

I am working with ArcGIS desktop 10.2

One of our major data supplier started data delivery using XML format using FTP. Packets of data updates are coming every day through FTP and stored in SQL.

I am looking for best practice to import data into geodatabases.

I have already used FME which has many shortcomings and I am not satisfied as its a lengthy process to write geometries. Also Importing XML to ArcMAP endup with error (AdvanceToElement: Unable to find xml element (workspaceDefinition not found in xml doc))


texts: What coordinate system is this data is in?

texts:
I have a CSV with latitudes, longitudes, and place descriptions, I'm trying to figure out the CRS so I can reproject the data.

The data has no accessible documentation, and appears to be dumped out of an Oracle database.

Plotting the coordinates on a cartesian grid shows all of the United states and some places outside of it.

What is the best way to figure the CRS out?

Here are a few examples from (juding by the FIPS codes) Brooklyn, Boston, and LA

Longitude, Latitude, State_Code, County_Code, Station_Location
73948800, 40725490, 36, 47, BTWN NORMAN AVE & DRIGGS AVE
71023639, 42424997, 25, 25, NORTHEAST EXPRESSWAY
118442961, 34042668, 6, 37, N/O OLYMPIC BL; S/O JCT RTE 2


texts:
I started using QAD plugin in QGIS, and try to create a polyline, based on the image of the Openlayers background map. After I finish digitizing the polyline and try to zoom in or out, it changes its actual position. Both the new layer and the map have the same spatial coordinates.
 The first picture is from the digitizing, and is in the right place
 And this one is after I zoom-out the image


texts: Using Google Maps, layers of polygons and calculate its intersections

texts:
I am developing a web application that uses Google Maps and it has polygons (circles and irregular polygons) drawed inside this map.

This polygons are stored into a DataBase (MySQL 5.7) as Polygons represented by the library JTS.

In the web application you can draw a polygon and then, a request is done to a server to check whether ths drawn polygon intersects or not with any other polygon stored in DB.

This process is quite slow and another approach is checking the intersection only in frontend. These is some solutions provided by other mates:

Using ArcGIS JS. I don't know if in this approach you could use Google Maps or another map is needed.
Using OpenLayers or Leaflet and maybe ArcGIS JS?
Using JSTS and Google Maps, or maybe, as I mentioned before, OpenLayers or Leaflet?

The problem I have using JTS is that represents circles (specifying radius and center point) is not the same as in Google Maps. So I'm not sure if using JSTS is a good idea.


texts: Best way to filter points given a starting point and max travel time

texts:
I have a database of about 1M geolocated addresses and I'm looking for a way to sort/filter these records based on a starting point and max travel time.

Ex:

"find every baker shop within 30 min of {where I am}"
"find every {whatever} within 60 min of {my work address}"
Optional: "find every {whatever} within 60 min of {where I am} AND 30 min of {my work address}"

Is there a known algorithm for this type of queries ?

It looks like PostGIS + pgRouting would fit the job, but I don't know where to start.

I'm looking for a server-side solution that I'll query with a PHP application.


texts: I have contours lines along with their elevations in ArcGIS. I want to create a polygon representation from those contour lines. How can i create polygons from those contours which may help me to represent them on my map.


texts: Arc GIS 10.1 - How to sort Layers ascending at TOC in "List by drawing order"

texts:
I have some 100 features in my ArcMap Project. I choose at TOC the View mode "List by drawing order".
How can I sort this list ascending or by name alphbetical?
Under the options field, is nothing like this.
The drawing order should be ascending by name.
Or how can I sort them on other way?


texts: Can i export BUS routes from OpenStreetMap as shapefiles to work on ArcMap?

texts:
The OSM has all the routing for the buses in my city, even showing the ID of the bus if you zoom close enough. I know how to download the shapefiles with the roads, rails, etc., from Geofabrik, but the bus routing isn't included there.

If I can't get access to this routes in ArcMap, I'll have to draw them one by one, which doesn't make that much sense since that work has already been done. Please, can someone help me?


texts: Problems using v.split in SEXTANTE toolbox to split lines

texts:
I got this error：

Or don't click "open output file after running algorithm", but there is no new shapefile created on my desktop.

How to solve this problem and split lines at a specified distance？

PS. I use QGIS 1.80 under WIN7-64 bits


texts: Printing label and associated pixel value from raster?

texts:
I have created a classified raster using the ArcGIS Pro image classification wizard. Unfortunately, I am unable to get the table of pixel intensities & associated labels to export to .csv.

When I export the raster, I get a summarized table of the classes and counts. However, I want to train a custom neural network with these images and want to compare the output from the classification raster (R G B individual values across the image).

Is there a way to do this?

I do not have access to the ESRI Imagery Analyst license so am unable to export the training for deep learning.


texts: Proper way to calculate a selection of GRIB2 data points?

texts: I have a GRIB2 file that I have successfully parsed and got an array of temperatures for a specific latitude and longitude. What would be the proper way to convert this array to a single value? I wasn't sure if taking the average might be the correct way to do this, so I wanted to double-check.


texts: Parsing Geoprocessing Service messages using ArcGIS API for JavaScript?

texts:
I am working on a JavaScript (JS) app where the user draws a polygon and then a python-based Geoprocessing (GP) Service task queries the parcels and generates mailing labels.  The GP task works great, but once the PDF is created I need to open the hyperlink in a new tab.

What I really need to get at is the return of my Python function:

# as a trivial example:
def myFunc(*args):
    # do something
    return output_url  # this is what I want to access!

So in my gp task, the function is returning the hyperlink url needed to download the report.  I do not think there is an out of the box way to get the gp result.  So the other thing I am thinking is to parse the messages returned from the GP Service.  I have looked through the help docs and do not see a sample how to get at the messages.  I'm new the JS development, so it's probably there and I'm just missing it.

My call to the gp task:

var gp = new Geoprocessor("http://myServer.com/arcgis/rest/services/GP_TASKS/MailingLabels/GPServer/MailingLabels");

function OpenInNewTab(url) {
    var win = window.open(url, '_blank');
    win.focus();
  };

var gpParams = {
    "Boundary": fs,
    "Site": "Paynesville",
    "Mailing_Type": "taxpayer"
  };
gp.execute(gpParams);
var URL = ""; #need to get this from result OR parse messages

gp.on("execute-complete", OpenInNewTab(URL));

And here is an example of the JSON from the GP service messages:

{
 "results": [

 ],
 "messages": [
  {
   "type": "esriJobMessageTypeInformative",
   "description": "Executing (MailingLabels): MailingLabels \"Feature Set\" Paynesville taxpayer"
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "Start Time: Wed Mar 25 11:29:45 2015"
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "Executing (MailingLabels): MailingLabels \"Feature Set\" Paynesville taxpayer"
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "Start Time: Wed Mar 25 11:29:45 2015"
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "Running script MailingLabels..."
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "Creating labels in Avery 5160 format"
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "\nCreated \\\\arcserver2\\wwwroot\\TempFiles\\MailingLabels_taxpayer20150325112946.pdf\n"
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "http://myServer.com/TempFiles/MailingLabels_taxpayer20150325112946.pdf"
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "Completed script MailingLabels..."
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "Succeeded at Wed Mar 25 11:29:47 2015 (Elapsed Time: 1.85 seconds)"
  },
  {
   "type": "esriJobMessageTypeInformative",
   "description": "Succeeded at Wed Mar 25 11:29:47 2015 (Elapsed Time: 1.85 seconds)"
  }
 ]
}

There is nothing in the "results", so since I am actually printing out the URL needed in the messages, I'm thinking I can parse these somehow.  So my questions are:

Is there a way I can get the return URL from my Python function back as a result from the GP Service?
If the answer to number 1 is no, how do I get at the result messages and step through them?  This is a Synchronous task so I am using the execute() method.


texts:
I would like to run this code that makes calculation at every loop between different input rasters and Itot that is a raster that in the loop remains always the same.
The code doesn't give me any error but it doesn't work because at the end the writing:'invalid layer gdal provider cannot open the gdal database' appear.
I guess the problem could be in how I treat Itot in the calculation, can someone help me?
Here is the code:
import qgis
import gdal
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

p=[0.84, 0.83, 1.03, 1.11, 1.24, 1.25, 1.27, 1.18, 1.04, 0.96, 0.83, 0.81]
Itot = QgsRasterLayer('/Users/macbook/Desktop/TESI/CAPITANATA PROGETTO/Itot.tif')

suffix_output=['th_jan', 'th_feb', 'th_mar', 'th_apr', 'th_may', 'th_jun', 'th_jul', 'th_aug', 'th_sep', 'th_oct', 'th_nov' ,'th_dec']
suffix_input=['Tjan', 'Tfeb', 'Tmar', 'Tapr', 'Tmay', 'Tjun', 'Tjul', 'Taug', 'Tsep', 'Toct', 'Tnov', 'Tdec']

outputpath = r'/Users/macbook/Desktop/TESI/PROGETTO/PET_thornthwaite/'
inputpath = r'/Users/macbook/Desktop/TESI/PROGETTO/raster_T_average/'

list=[0,1,2,3,4,5,6,7,8,9,10,11]
for i in list:
  inputrasterfile = QgsRasterLayer(inputpath + suffix_input[i] + ".tif")

  entries = []
  ras = QgsRasterCalculatorEntry()
  ras.ref = 'ras@1'
  ras.raster = inputrasterfile
  ras.bandNumber = 1
  entries.append( ras )

  new_path = outputpath + suffix_output[i]+ '.tif'

  th = str(p[i]) +'* 16* (10* ras@1 / Itot)^(0.5 + 0.016 * Itot)'

  Pet = QgsRasterCalculator(th, new_path, 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries )

  Pet.processCalculation()
  iface.addRasterLayer(new_path)


texts: Extracting features to feature classes from a REST service

texts: I've seen this done before but I can't figure that out how to do it. How do you extract features from a REST service to a feature class (or shapefile)? I believe it involves adding a REST service link to a map in ArcGIS Online (AGOL) then opening in ArcGIS Desktop. Anybody know how to do this?


texts: General steps to create watershed boundary in ArcGIS Spatial Analyst?

texts:
I am an ArcGIS beginner.

I have a CD of GIS files of a city and intend to work on watershed management.

There are different file names in different formats on CD. Like .shp, .shx, .dbf, .sbx. and I have no idea which one to choose to start working on.

As there is a file named contour line (in all mentioned formats) I would guess I need to start with that one. I am not sure, though. I also found out that I should use Spatial Analyst Extension but in which step I do not know.

What are the main steps, and which file to start with to end up with a 3D file of watershed and its boundaries?


texts: Exporting WKB (geometry) to .txt using ArcPy?

texts:
How can I export WKB (geometry) to text file using arcpy (ArcGIS)?

When I use the examples that I found on the Internet, it doesn't work. Below it's the code I'd wrote:

## Import arcpy module

import sys
import string
import os
import arcpy

## Global variables:

mxd = arcpy.mapping.MapDocument("current")

IN_FC = arcpy.GetParameterAsText(0)
OUT_Addre = arcpy.GetParameterAsText(1)
OUT_Name = arcpy.GetParameterAsText(2)

## Set progressor
lines = arcpy.GetCount_management(IN_FC)
number = int(fc.getOutput(0))

arcpy.SetProgressor("default","Exporting WKB geometry",0,number,1)

## for row in input
with arcpy.da.SearchCursor(IN_FC, ["OID@", "SHAPE@WKB"]) as cursor:
    for row in cursor:
## First feature
        elem = IN_FC.getFeatures().next()
        geom = elem.geometry()

## Export into wkb format
        wkb = geom.asWkb()

## Import to ogr geometry
        from osgeo import ogr
        geom_ogr = ogr.CreateGeometryFromWkb(wkb)

## Import to Shapely geometry

        from shapely.wkb import loads
        line = loads(wkb)
        ls = wkb.encode('hex')

## output

        output_file = open(OUT_Addre & OUT_Name & ".txt", 'w')
        output_file.write('\n' & ls.encode('utf-8'))
        output_file.close()

Some errors message:

  import ogr or wkb = geom.asWkb() "doesn't exist"


texts:
Can you please help me out in Installing pgRouting in Windows 7 64-bit version.

I have been trying it for days. Since I am completely new to this please provide me the initial requirements and steps in installing pgRouting.


texts: I'm trying to create a timeseries of raster data (one layer per year) by polygon (Administrative regions). Is there any way in QGIS to perform Zonal Statistics for multiple years (without having to perform the zonal statistics function each time)? Or something along the lines of the Point Sampling tool but for polygons?


texts: Nearest neighbor distance in km over contiguous US

texts:
I have two sets of line features in QGIS (rivers & roads - I'd like to find the road segment nearest to each river segment & the distance between them in km) across the contiguous United States. I've tried using the nearest neighbor join plugin (NNJoin) with both WGS 84 and the USA contiguous equidistant conic (ESRI: 102005) projections, but I get distances between 0 & 50, which I'm assuming is in degrees. Should I try a different projection or is there a way to manipulate the NNJoin plugin?
Edit: I have also tried keeping the river segments in WGS 84 & only putting the road segments (join layer) in a projected CRS (as per the help page on the plugin). However, this reprojects the road layer to a totally different location from the river feature & also gives clearly incorrect distance values.


texts: layer styling switches after zooming in with no reason

texts:
Does anybody have an explanation how it is possible?

I have a strange trouble on the following map. The surrounding land layer, which is only a one polygon layer and is styled normally like:

var stylePoly = new OpenLayers.StyleMap({
        fillColor: 'white',
        fillOpacity: 0.85,
        strokeColor: 'black',
        strokeWidth: 1
        });

is switching to the non-surrounding area when zooming to a certain level. It only happens when the focus is on a bbox where the surrounding land is to be seen.

What should I change to get the appearance constant? I do not really understand this behavior. Maybe somebody had the same problem already...

The data is gml and the visibility of the layer is constantly "true"...

Layer is created like:

var protocolUmland = new OpenLayers.Protocol.HTTP({
        url: 'http://bildungsatlas.map-dienst.de/gml/FHH_Umland.gml',
        format: gml
        });

var Umland = new OpenLayers.Layer.Vector('Umland',
    {
    strategies: [new OpenLayers.Strategy.Fixed()],
    protocol: protocolUmland,
    styleMap: stylePoly,
    visibility: true
    }
);


texts: Script Error 'Invalid Pointer' in Spatial Analyst Tool

texts:
I have ArcGis 10.3 (previously 10.2 version) running on Windows 10 (previously Windows 8). My GIS licence is provide by my company where I have an intership, so I can't just uninstall/install back the logiciel.

So, I have this script error which appears often when I want to use Spatial Analyst Tool or Geoprocessing Tool:

  Line 51, Char 4, Error Invalid Pointer, Code 0, URL
  file:///C:/Users/username/AppData/Roaming/ESRI/Desktop10.3/ArcToolbox/Dlg/MdDlgContent.htm

I look at tons of articles online and did all their proposed solutions (put Internet Explorer as default browser, enable Active X control in security, clear roaming and local history of ArcTool Box in my AppData, upload a 'new patch version' of MdDlgContent.htm in my ArcGIS file provide by technical support from ESRI) and NOTHING IS WORKING.


texts: Randomly distributed cases according to census data

texts:
I'd like to randomly place a point within a geographic boundary with the caveat that I would like to use census data so that points are randomly placed according to population density (i.e. more likely to be placed where there is higher population density).  I would ideally like to do this using the R programming language.

I am not a GIS person, so I am not entirely certain of how to do this. I have some census data and  I have my boundary co-ordinates (as a shapefile, and this I know how to use), but how I go about combining the two and more so, placing a point randomly according to census data is where I get a little stumped.


texts: Transfering real time data from SQL Server to ArcMap without using PI Integrator

texts:
I am new to GIS. Right now I am working on my senior project in which I have to find the fault in transmission line and show its location on a map using ArcMap 10.2. I have some readings of voltage and current in Microsoft SQL Server 2016. I want to transfer these values in REAL TIME to ArcMap 10.2.
Is it possible to do it directly, or do I have to use PI Integrator for Esri ArcGIS?


texts:
I can not get text element to update?

Text ELement has been added "Subdivision" Still not updating?..

class ButtonClass2(object):
"""Implementation for FemaMapper_addin.button (Button)"""
def __init__(self):
    self.enabled = True
    self.checked = False
def onClick(self):
    #set the map document reference
    mxd = arcpy.mapping.MapDocument("CURRENT")
    #set the dataframe reference using the map document and the first data frame in the list
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    #set the layer reference using the map document reference and the second layer in the list of layers
    layer = arcpy.mapping.ListLayers(mxd)[1]
    #heading of the map
    map_title = "Parcel Map"
    #create a search cursor to get relevant information about a parcel
    sc = arcpy.SearchCursor(layer)
    #for loop to get the selected data
    for row in sc:
        #select the relevant fields to be put on the map
        Sub = row.getValue("SUBDIVISIO")

    #text element of the map
    elm = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "Heading")[0]
    elmSub = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "Subdivision")[0]

    #font of the heading
    elm.text = '<FNT name ="Arial" size="32">' +  map_title + '</FNT>'
    elmSub.text = '<FNT name ="Arial" size="14">' + 'Subdivision: ' + Sub + '</FNT>'

    #Change the view to layout view
    arcpy.mapping.MapDocument("CURRENT").activeView = "PAGE_LAYOUT"
    #refresht the active view
    arcpy.RefreshActiveView()
    #delete search cursor and row variables
    del sc
    del row


texts: convert GPS coordinates to SRID coordinates in spatialite

texts:
The following query fails for me:

select account_nu, srid(geometry), AsText(Geometry)
from tbl
where contains(geometry,MakePoint('35.16604','-101.8635',4326));

I would like to find shapes that contain a specific point.  All of my shapes are stored under srid 6678.  I assumed the GPS coordinates are stored as 4326, but that's not the case.

How do I convert GPS coordinates to coordinates that I can use in spatial lite?


texts: QGIS Clip function in Raster layer Style Single band pseudocolour not working?

texts:
I've been using  QGIS for a while and recently upgraded to 2.4, but now find that displaying raster layers is difficult. In particular, the clip function (tick box bottom left of "Style">"Singleband pseudocolor" menu) does not seem to be working. This makes it difficult to display a number of layers together.

To specify the problem: I wish to display cells with values ranging between 20-200 on a layer where values range between 1 and 200, and would have added values 20 and 200 into the window and ticked the "Clip" box in previous versions, but this no longer works. All cells with values 1-200 are displayed.


texts: How does gridmetrics handle IF condition in user defined functions?

texts:
I noticed a kind of "inconsistency" in the values returned from two user-defined functions. I tested on a tile of 200mx200m to return a raster of pixel size 50m (16 pixels)
In the following function, if the condition is satisfied, the function returns ONLY one value.
funcky <- function(x)
{
  if (max(x) < 22)
     return(max(x))
}
grd <- grid_metrics(lasca, funcky(Z), 50)

When the condition is not satisfied, NA values are returned.
> as.array(grd[[1]])
, , 1

      [,1]  [,2]  [,3] [,4]
[1,] 19.92 21.90    NA   NA
[2,]    NA 21.46 21.89   NA
[3,]    NA 21.73    NA   NA
[4,]    NA    NA    NA   NA

I expected the same logic to be applied when the function returns a list of values. For example,
funcky <- function(x)
{
  if (max(x) < 22)
    return(list(max(x), mean(x), min(x)))
}
grd <- grid_metrics(lasca, funcky(Z), 50)

But the values are not NA, they are 0s as shown below.
> as.array(grd[[1]])
, , 1

      [,1]  [,2]  [,3] [,4]
[1,] 19.92 21.90  0.00    0
[2,]  0.00 21.46 21.89    0
[3,]  0.00 21.73  0.00    0
[4,]  0.00  0.00  0.00    0

, , 2

         [,1]     [,2]     [,3] [,4]
[1,] 6.555234 7.434102 0.000000    0
[2,] 0.000000 5.038366 5.733465    0
[3,] 0.000000 8.699386 0.000000    0
[4,] 0.000000 0.000000 0.000000    0

, , 3

      [,1]  [,2]  [,3] [,4]
[1,] -0.83 -0.04  0.00    0
[2,]  0.00 -0.07 -0.09    0
[3,]  0.00 -0.07  0.00    0
[4,]  0.00  0.00  0.00    0

If the metrics can have negative values (skewness, for example), then these 0s could potentially be misleading, as is the case in the third layer. Is this behaviour alright?
Adding an else {return(NULL)} didn't make any difference. And adding else {return(list(NA,NA,NA))} generated an error:
Error: Column 1 of result for group 2 is type 'logical' but expecting type 'double'. Column types must be consistent for each group.


texts: OpenLayers weird problem with GoogleMaps base layer and WMS overlay

texts:
I'm trying to overlay a WMS layer over GoogleMaps base one. The WMS
layer is defined like this:

 var catastro_layer = new OpenLayers.Layer.WMS(
           "WMS Catastro",
           http://ovc.catastro.meh.es/Cartografia/WMS/ServidorWMS.aspx?SERVICE=WMS&VERSION=1.1.1,
           {
                   layers: "catastro",
                   format:"image/png",
                   srs: "EPSG:3785",
                   transparent: false
           }
   );

GoogleMaps layers are defined just like here
http://openlayers.org/blog/2010/07/10/google-maps-v3-for-openlayers/.

As result, my WMS layer appears as one base layer more. So, I can
select it instead of the other base layers, not over them. This
screenshots shows what I mean:
http://dl.dropbox.com/u/6599273/capa_wms1.png

The WMS requests are ok:

http://ovc.catastro.meh.es/Cartografia/WMS/ServidorWMS.aspx?SERVICE=WMS&VERSION=1.1.1&LAYERS=catastro&FORMAT=image%2Fpng&SRS=EPSG%3A3785&TRANSPARENT=false&REQUEST=GetMap&STYLES=&EXCEPTIONS=application%2Fvnd.ogc.se_inimage&BBOX=-1721973.3737,3228700.0733875,-1702405.4944625,3248267.952625&WIDTH=256&HEIGHT=256

But this isn't what I need. I want a GoogleMaps base layer and a
transparent WMS overlay.

Now, if I change the WMS layer defintion replacing transparent: false
by transparent: true:

 var catastro_layer = new OpenLayers.Layer.WMS(
           "WMS Catastro",
           http://ovc.catastro.meh.es/Cartografia/WMS/ServidorWMS.aspx?SERVICE=WMS&VERSION=1.1.1,
           {
                   layers: "catastro",
                   format:"image/png",
                   srs: "EPSG:3785",
                   transparent: true
           }
   );

My WMS layer appears as a overlay:
http://dl.dropbox.com/u/6599273/capa_wms2.png, just what I want.

But now, the requests are incorrect:

http://ovc.catastro.meh.es/Cartografia/WMS/ServidorWMS.aspx?SERVICE=WMS&VERSION=1.1.1&LAYERS=catastro&FORMAT=image%2Fpng&SRS=EPSG%3A900913&TRANSPARENT=true&REQUEST=GetMap&STYLES=&EXCEPTIONS=application%2Fvnd.ogc.se_inimage&BBOX=-1721973.3737,3228700.0733875,-1702405.4944625,3248267.952625&WIDTH=256&HEIGHT=256

As you can see, the epsg:3785 has been replaced by epsg:900913. But
the only change I've made is "transparent: false" by "transparent:
true". Why is this happening?

Many thanks in advance, and best regards


texts: How to use basic authentication for WMS layer in Leaflet

texts:
I use Leaflet.

A governmental agency provides a WMS layer with basic authentication. I first added and published it on my Geoserver and it worked, but I could not access the different styles. They suggested me to connect directly to their URL without going via Geoserver. How do I provide the credentials (username / password) for the authentication in a secure way when connecting directly to their GetCapabilities URL?

var wmsLayer = new L.TileLayer.WMS('http://maps-ver.lantmateriet.se/fastighet/wms/v1.1?request=GetCapabilities&version=1.1.1', {
                layers: 'granser',
                styles: 'morkbakgrund',
                format: 'image/png',
                transparent: true
            }).addTo(mymap);


texts: R and ArcMap not showing the same raster values

texts:
I imported one of my rasters into ArcMap to see how it is displayed in this program and noticed that it's values do not correspond to the same values I had in R.

Here are some screen caps:

In R:

Info about the raster (.asc format) in R:

  Formal class 'RasterLayer' [package "raster"] with 12 slots
  ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
  .. .. ..@ name        : chr "D:\\Users\\cost_guinea_29N.asc"
  .. .. ..@ datanotation: chr "FLT4S"
  .. .. ..@ byteorder   : chr "little"
  .. .. ..@ nodatavalue : num -Inf
  .. .. ..@ NAchanged   : logi FALSE
  .. .. ..@ nbands      : int 1
  .. .. ..@ bandorder   : chr "BIL"
  .. .. ..@ offset      : int 0
  .. .. ..@ toptobottom : logi TRUE
  .. .. ..@ blockrows   : int 1
  .. .. ..@ blockcols   : int 909
  .. .. ..@ driver      : chr "gdal"
  .. .. ..@ open        : logi FALSE
  ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
  .. .. ..@ values    : logi(0)
  .. .. ..@ offset    : num 0
  .. .. ..@ gain      : num 1
  .. .. ..@ inmemory  : logi FALSE
  .. .. ..@ fromdisk  : logi TRUE
  .. .. ..@ isfactor  : logi FALSE
  .. .. ..@ attributes: list()
  .. .. ..@ haveminmax: logi FALSE
  .. .. ..@ min       : num Inf
  .. .. ..@ max       : num -Inf
  .. .. ..@ band      : int 1
  .. .. ..@ unit      : chr ""
  .. .. ..@ names     : chr "cost_guinea_29N"
  ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
  .. .. ..@ type      : chr(0)
  .. .. ..@ values    : logi(0)
  .. .. ..@ color     : logi(0)
  .. .. ..@ names     : logi(0)
  .. .. ..@ colortable: logi(0)
  ..@ title   : chr(0)
  ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
  .. .. ..@ xmin: num -176511
  .. .. ..@ xmax: num 654315
  .. .. ..@ ymin: num 790241
  .. .. ..@ ymax: num 1408105
  ..@ rotated : logi FALSE
  ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
  .. .. ..@ geotrans: num(0)
  .. .. ..@ transfun:function ()
  ..@ ncols   : int 909
  ..@ nrows   : int 676
  ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA
  ..@ history : list()
  ..@ z       : list()

In ArcMap:

Info about raster in ArcMap:

Same layer, same source, same format, but different values.. why?
R seems to display the correct values whereas ArcMap doesn't.


texts: Exporting PDF from QGIS with print composer without destroying fonts

texts:
So I'm making this nice map with a lot of labels which are going to end up in a book. Since the labels doesn't end up exactly where I want them I'd like to move them some in Illustrator after exporting them using the Print composer. The only problem is that when exporting to PDF the text is not text anymore, QGIS transforms it into paths instead of leaving the font alone.

Does anyone know a workaround so the labels stays as fonts when opening the image in Illustrator?


texts: Point shapefile missing spatial reference information

texts:
I am using ArcGIS 10.2.2. I've got a shapefile with points representing spot elevations in my study area. This file has been provided by a third party institution.

The coordinate system is missing spatial reference, as ArcGIS informs me as soon as I feed the file into ArcMap. I tried to use the Define Projection tool, but nothing happens, i.e. when I add other feature classes created by me (and that matches one another) the shapefile in question does not 'match' the rest of the layers.

As I tried to explore the extents under source tab in properties, the coordinates are shown as:

left 3199745.231567 ??
right 3219058.255825 ??
top 1064279.877094 ?? and
bottom 1041353.944159 ??

When I contacted the source party, they informed me that actual coordinate system of the file generated is UTM Zone 43 N but as file was originally generated in CAD and then exported to GIS, its .prj file is missing.

Moreover coordinates shown in the problematic shapefile are not corresponding to actual coordinates in the right coordinate system (UTM Zone 43 N). For example: 'true' coordinates of a point 319907.40 E / 3720430.82 N, corresponding coordinates in the problematic  3212077.85 E / 1048297.355 N; so, the first is missing the one digit six digit as compared to seven, the second is showing different coordinates all together.

I would like to know if there is any way to fix this problem?

As I received that data, file with extension .prj is missing in the folder.


texts:
How can I load .SID files in ArcEngine?   I think its possible, but I can't find any links.

With Kirk's help I got it here's the code:

  public void add_SID_file(string dir)
        {
            RasterWorkspaceFactory raster = new RasterWorkspaceFactoryClass();
            IRasterWorkspace rasterWorkspace = raster.OpenFromFile(dir, 0) as IRasterWorkspace;

            if (rasterWorkspace != null)
            {
                IWorkspace w = rasterWorkspace as IWorkspace;
                IEnumDatasetName dsNames = w.get_DatasetNames(esriDatasetType.esriDTRasterDataset);
                IDatasetName datasetName = dsNames.Next();
                List<string> names = new List<string>();
                while (datasetName != null)
                {
                    names.Add(datasetName.Name);
                    datasetName = dsNames.Next();
                }
                if (names.Count > 0)
                {
                    IRasterDataset pRasterDataset = rasterWorkspace.OpenRasterDataset(names[0]);
                    IRasterLayer pRasterLayer = new RasterLayerClass();
                    pRasterLayer.CreateFromDataset(pRasterDataset);
                    mapControl.Map.AddLayer(pRasterLayer);
                }
            }


texts: Resolving Geoserver ERROR: Operation on two geometries with different SRIDs?

texts:
I am working on GeoServer Version 2.4.2 . I have one layer 'beat' in PostgreSQL database. I have publish layer in geoserver but its not opening OpenLayers layer preview link. It gives error as

Internal error occurred java.io.IOException null ERROR: Operation on two geometries with different SRIDs

I have checked all previous question but i didn't understand how to resolve this. I try exporting the layer and it export successfully in shp and render in QGis. Also it export in geojson succesfully so why its not rendering in layer preview i am confused. I tried following queries but it shows proper results as

select * from geometry_columns where f_table_name = 'beat';
    f_table_catalog  f_table_schema  f_table_name  f_geometry_column  coord_dimension  srid   type
    ---------------  --------------  ------------  -----------------  ---------------  -----  ------------
                     nmmccommon      beat          the_geom           2                32643  MULTIPOLYGON

select distinct(st_srid(the_geom)) from beat;
st_srid
-------
32463

Distinct srid is also single then at what point Operation on two geometries with different SRIDs  error occure. Geosever can export it shp also successfully render in QGis. I know its repitative question but i didn't found a useful solution which suits in my case.


texts:
I asked a Map related question on Stack Overflow:  https://stackoverflow.com/questions/9948710/offline-ovi-map-in-desktop-application

To provide more details, we are focussed on developing a commercial desktop application that will not have any internet connectivity. Though open source OpenStreetMaps seems like a good option, a commercial solution is preferred.


texts: I'm using ArcMap 10.1. I want to create buffers that represent CO2 equivalent from landfills.  I used graduated symbols to label each landfill site, based on amount on incoming waste per year.  Will the size of the symbols impact the size of my buffers? Or do buffers always measure diameter from the centroid of the point?


texts:
I'm a newbie to GIS and QGIS (please be nice to me).  I'm working on a map for a garden/arboretum catalogue using QGIS 2.6.1.  The arboretum has been surveyed using GPS, and the specimen trees and plants become a point layer labelled with a reference number, which points to a description in the main catalogue.

Because of the site layout and orientation, the map needs to be rotated to fit comfortably on a rectangular canvas (paper, screen).  Fine, I can do that in print composer.  But that rotates the labels too, and for maximum legibility I want them to stay horizontal.

I've been experimenting with label placement. "Around point" gives significantly better results than "Offset from point".  Using this method, the only way I've found to apply a compensating rotation to the labels (so that when the map is rotated they end up back horizontal) is to add a "rotation" column to the attribute table and use data defined rotation.  But this to some extent messes up the label placement.  So my first question is:  is it possible to get the placement engine to run (again?) after the labels have been rotated?

If not, I would like to be able to reposition a few of the labels manually, but leave most of them automatically placed.  But if I add coordinate x and y columns to the attribute table and try to use data defined coordinates, then label rotation stops working if there is a NULL value in either or both of the x and y coordinate columns.  I can't leave most of the points to be automatically placed and just fix a few - effectively I have to do manual placement for all of them if I want to do it for any of them.  Since there are approx. 200 of the things, that's a major chore.  Have I missed something?

Some of this might be easier if the "Move label" button would work.  But however I've fiddled - with editing on/off, data defined placement on/off and so on - I've never seen the thing not greyed out.  What do I need to do to get this working?

The best I've been able to do is have a data defined distance override for some points.  This is rather hit and miss, because there's not much control over where the label will actually go, and it has unpredictable knock-on, domino effects on other labels some of which may end up worse than they were before.  But it doesn't affect label rotation.


texts: How can I find points inside a polygon in ArcGIS server?

texts: I have a feature class of points. My client wants to draw a polygon on the map to select the subset of points inside the drawn polygon. Is there any tools in ArcGIS API, which I can use it for this purpose? if there is a sample for ArcGIS API for Flex, that would be great.


texts: Adding spatial information to a table, using spatial join but getting null columns

texts:
I know this could be a duplicate (Merging attributes from point with containing polygon?), but the answer in the link is not clear. I have a very similar task where I want to add district information in a shapefile to a table with attributes of schools in Shanghai area.

This table has information on more than 2000 schools and XY coordinates, but no information on district. I want the polygon information of the district to be added to the school attributes table, so that I can identify the districts each school is located.

As the link directs, I tried spatial join, but the information of the shapefile were not added appropriately to the school attributes table. As in the picture, all shapefile information was added as .

I set the target feature as school attributes table with XY coordinates (sh_np2_Layer) and the join features as the district polygon shapefile (shang_dis_merged), and join operation as join one to many, because I want to join district information to each school point.

I'm bit lost on the parameters part, how should I set it correctly so that I can join the layer and the shapefile correctly? I am especially lost on setting the 'output field' , 'merge source', and 'match option'.


texts:
I have a folder full of about 2000 files and I would like to generate a kind of key map for each file so people don't need to just randomly load each one in order to find what is what.

I have the following python code using ogr but OGREnvelope isn't really a OGRGeometry so I can't add it to a shape file or export it as wkt to import in QGIS.

for filename in glob.glob(path + "*.shp"):
   ds = ogr.Open(filename)
   layer1 = ds.GetLayer(0)
   #print filename + " : " + str(layer1.GetExtent())
   extent = layer1.GetExtent()
   # I need to be able to export the extent region as a polygon.

I have been looking around but just can't seem to find anything.  I really thought something like this would be easy.


texts: Opening marine phytoplankton NPP data in QGIS

texts:
How can I open these file in QGIS? I am beginner at GIS softwares. My aim is to obtain marine phytoplankton NPP data.

These data are available at oceancolor.gsfc.nasa.gov and www.nodc.noaa.gov, but QGIS won't open the files as simply as I thought.


texts: Google map Integration with openlayers 3 features displacement error on map zoom and drag

texts:
If I add google map as in this link (http://openlayers.org/en/v3.0.0/examples/google-map.html), when I zoom or drag the map, features also displacing from their position and coming back to its actual position. How can I add Google Map as base map without displacement of features?

I don't want to use: https://github.com/mapgears/ol3-google-maps, because of this library I am getting many errors in the application and it doesn't support all functionalities of openlayers 3.


texts:
How to rank polylines according to their similarity to a ref polyline?

Update 1: simulations show that simple $\sqrt{\sum{(case-ref)^2}}$ can find better similarity between polylines (as set of points) than Hausdorff distance! but not satisfied yet! (right-click >> view image >> for full resolution)


texts: Manually setting spatial extent of custom map in QGIS 3.14

texts:
Due to customization wishes I created a map in Adobe Illustrator which means my .tiff file has no metadata like CRS, extent etc. when I export it.
To manually set the spatial extent I would need an option like it used to be in the properties menu (How to edit a vector layer extent?) but it was removed.
I can't find a similar option in a newer QGIS version (3.14 or higher) for raster data. Though the menu to set the extent metadata is still there it is only for documentation purpose and does not actually change the size of the map on the canvas. If you know how to achieve this while keeping RGB values, please help me.


texts: Is there any way to mass rename or import field names in Labels in Symbology?

texts:
I need to load an .avl that sets the symbology for me, but I also need to give the unique values their name, which is not what the Symbology->Category->Value Field uses, but another field in that layer (or in another, but that's just a join...) I do not know if there is a way to import other fields for labeling.

I am using 10.2.2 version of ArcMap.

I only want to change the labels, not one by one, but to change it to another field. The coloumn marked with the arrow should go like "a,b,c,d...".


texts: Does the performance of Postgis function ST_DWithin degrade with larger distances?

texts:
I am using the Postgis ST_DWithin function like this:

SELECT *,
    ST_X(geometry) as longitude,
    ST_Y(geometry) as latitude,
    ST_Distance_Sphere(
        ST_SetSRID(
            ST_MakePoint(-118.478255,34.303046999999999),
            4326),
        geometry) / 1609.34 as DISTANCE
FROM  poi
WHERE ( category =10980400 )
    AND ST_DWithin(geometry::geography,
        (ST_SetSRID(
            ST_MakePoint(-118.478255,34.303046999999999),
            4326)::geography
        ),
        80467.2,
        false)
ORDER BY DISTANCE
LIMIT 10

I am not sure if the performance degrades with increasing distance like 5, 10, 20, 40 miles. Please help me understand if it degrades then what can be the reason ? Is it going to be exponential or will it depend on the number of points ?

I am using it for finding nearest points (like ATMs, restaurants) within a radius.

I have indexes available for geometry and geography. When I "explain analyze" the query I can see the usage of the spatial index. I am using geography because I am getting user input in meters, kilometers etc.

Is there a way by which I can stop degradation over large distances ? We are supporting maximum 50 miles.

The problem with using KNN is that it degrades in those cases where it doesn't find 10 results and keep on searching till any distance. Initially we were doing a subquery using KNN and then filtering the results on distance outside the subquery. It timed out in those cases where it started searching for let's say an ATM and kept on searching till it found 10. If the radius by user is only 1 KM then also it kept searching outside that radius till it found 10.


texts:
Given an area id, is there a way to get size of this area from Overpass API?

For instance, I would like to know how many square kilometers New York has.


texts:
I am trying to follow a tutorial from here http://qgistips.spatialthoughts.com/2012/12/find-neighbor-polygons-in-layer-in-qgis.html#more

I have downloaded the relevant python script and edited it as required. When I try to run the script using the Python console in QGIS >Plugins>Python Console I get the following error message.

ImportError: No module named shapely.wkb
Start of relevant script

from PyQt4.QtCore import *
from shapely.wkb import loads
# Replace the below value with the field containing name or id of the feature
# For example, if your field is called name then change the line below to
# name_filed = 'name'

The module shapely.wkb is called on the second line of the script so the failure appears to be occurring as the second line of code is executed.

I have searched for a Python module with that name in Plugins>Fetch Python Modules
Is this something that needs to be installed separately? From what I have read it should be installed with an installation of PostGis which I though got loaded when I used StackBuilder to install PostGreSQL.

Would installing installing Python 2.7.3 on my MAC Book Pro help?


texts: QGIS Modeler: Multiple layers going through the same process

texts:
Current I am trying to build a process model that edits fields and attributes of a collection of vector layer, the thing is though I can only select one vector layer and time to go through the "Refactor Fields" Algorithm. But because all three of the Vector layers have the same kind of attribute table and attributes. I want all of these three layers, going through the process at once. Not that we have to do the process three times but only once. How do we do this?

Note: this model is still a build in process, it is the refactor fields algorithm if you cannot understand dutch.


texts: ArcPy Editor commits in batches of a thousand records dropping some records

texts:
I am copying data from one database to another using arcpy cursors.  This is an Oracle enterprise geodatabase against a versioned dataset.  Using Python 2.7 and ArcGIS 10.1.

The following code has worked in the past:

def copy_simple(target_workspace, source_fc, target_fc):
    dsc = arcpy.Describe(source_fc)
    fields = dsc.fields
    fieldnames = [field.name for field in fields if
                  field.name != dsc.OIDFieldName
                  and field.name != 'GLOBALID'
                  and field.name != 'SHAPE'
                  ]
    fieldnames = fieldnames + ["SHAPE@"]

    edit_target = arcpy.da.Editor(target_workspace)
    edit_target.startEditing(False, True)
    edit_target.startOperation()
    with arcpy.da.SearchCursor(source_fc, fieldnames) as sCur:
        with arcpy.da.InsertCursor(target_fc, fieldnames) as iCur:
            for row in sCur:
                inserted_oid = iCur.insertRow(row)
    edit_target.stopOperation()
    edit_target.stopEditing(True)
    del edit_target

However, while there are 17,592 records in the source feature class, there are only 17,000 records in the target at completion.  It seems that the editor is batching the commits in 1,000 increments and the last 592 are being dropped.

I found this post on GeoNet, but using the with did not fix the problem.

This code has worked in the past without issue for different databases and feature classes.  I am not sure what environmental setting would cause this odd behavior.

Update:

The following code that uses an append approach does capture all the records:

def append_simple(target_workspace, source_fc, target_fc):
    edit_target = arcpy.da.Editor(target_workspace)
    edit_target.startEditing(False, True)
    edit_target.startOperation()
    arcpy.Append_management(source_fc, target_fc, 'NO_TEST')
    edit_target.stopOperation()
    edit_target.stopEditing(True)

This might have to be the way I go, but it would be great to try and understand what is the cause of the cursor not being complete.


texts:
In C# how can i convert a Lat/Lon coordinate to Vicgrid94 which is described http://spatialreference.org/ref/epsg/3111/ and also http://www.dse.vic.gov.au/__data/assets/pdf_file/0012/117102/VICGRID94MapProjectionInformation.pdf. The second link gives math examples of how to convert but is a little to much for me.

I am aware of JavaScript library's and the like, but i would prefer to do this in C#


texts: Intersect two layers and get resulting features

texts:
I have X layers on a map and I need to intersect chosen layers (two per time) and color resulting features.

I'm trying to get this code working

// get first feature (index 0)
ESRI.ArcGIS.ADF.Web.Geometry.Geometry adfFeature =
                m_firstLayer.GeometryFromRow(m_firstLayer.Rows[0])
as ESRI.ArcGIS.ADF.Web.Geometry.Geometry;

// THE FOLLOWING LINE RETURNS NULL
ESRI.ArcGIS.Geometry.IGeometry featureInterface = adfFeature as ESRI.ArcGIS.Geometry.IGeometry;

ESRI.ArcGIS.Geometry.ITopologicalOperator topoOp = adfFeature as ESRI.ArcGIS.Geometry.ITopologicalOperator;

How to use the IGeometry interface with an ADF Geometry object?

I can't really find samples to intersect features between two layers, and it's a pity that Spatial Joins are just a arcgis desktop function, I surely could use them!

Help please


texts: Snap to crossing between lines in OpenLayers

texts:
Is it possible to let the snap interaction snap to the crossing of two different lines (or polygons) in OpenLayers?

It can snap to vertices and lines between vertices, but I can't see a way to snap to the crossing of different geometries.


texts:
When I run certain scripts (contains tkinter or similar python packages) ArcGIS crashes. When I load or copy/paste the same code and run it in a Python window it runs without error.

Is there a problem in ArcGIS? I am using 10.4.1.


texts: How to insert and transform geometry fields in PostGIS?

texts:
I have a source table with a geometry(Point,27700) column, and a target table with a geometry(Point,4326) column.

I am using the following Python code to join the two:

q = "SELECT ST_SetSRID(ST_AsText(geom), 4326) FROM source WHERE id=%s"
cursor.execute(q, (id)
geom = cursor.fetchone()
if geom:
    query = "INSERT INTO target(geom) VALUES(%s)"
    cursor.execute(query, (geom))

However, if I then extract the target table using the following command:

 ogr2ogr -f GeoJSON out.geojson "PG:host=localhost dbname=mydb" -sql "SELECT * FROM target"

The geom fields still look like 27700:

"geometry": { "type": "Point", "coordinates": [ 354151.0, 193500.0 ] } }

What am I doing wrong?


texts: I'm trying to get a Height Profile for a certain route using the profile tool (terrain profile). Does anyone know how i can get the slope in percent and in degrees in the table of the profile tool so that i can export it?


texts:
Attempting to import and convert Excel spreadsheets to tables then add XY Coordinates to the tables. The tool runs and provide the output tables below, but I cannot get the second part of the tool to run properly and populate XY coordinates. I know that I can manually add the XY coordinates, but having the tool automatically create them is necessary for the purpose of this tool.

When I run the tool, an error message say it cannot find the tables in the the specified file path.

Is it possible in the same tool to add fields to an EXCEL to TABLE conversion?


texts: Dismantling and reassembling polyline by joining each point to its nearest neighboring point using ArcGIS Desktop?

texts:
How can I reconfigure a PolylineMZ feature class in a way as shown in the image below?

I have a feature class (PolylineMZ) which contains x-sections. Some x-sections (polylines) are badly constructed (as shown below in figure "a"), so I want to dismantle and reassemble them again in the right order. By "right order", I mean that I want to join the points that are closest to each other (as shown below in figure "b").

Any idea for an existing tool? Preferably in Arc 10.


texts: pgr_drivingdistance returning no results

texts:
I've imported the latest North America OSM extract into postgres using osm2po so that I can use the pgr_drivingDistance function of pgRouting to create a drive-time polygon.  I'm at the final step, but the function is returning no results even though it definitely should.

I use a bounding box to narrow down the possible results inside the query consumed by `pgr_drivingdistance.  My first debug step was to see what that is returning.
My inner query (altered temporarily to include the geom):

SELECT rn.id AS id,
rn.source::int4 AS source,
rn.target::int4 AS target,
rn.cost::float8 as cost,
rn2.geom_way
FROM public.na_2po_4pgr rn
left join public.na_2po_4pgr rn2 on rn2.id = rn.id
where ST_Intersects(rn.geom_way, st_geomfromtext('POLYGON((-79.5670159782642 36.3965968758439,-79.5670159782642 36.6695646217025,-79.2319617622083 36.6695646217025,-79.2319617622083 36.3965968758439,-79.5670159782642 36.3965968758439))', 4326))

This query returns 10230 rows, and looks like:

That's correct, and is the correct region, and obviously, plenty of nodes.

I've got a start node, and picked it's ID, which I've verified is present in the returned rows:
14232131

So finally, I put it all together, and run pgr_drivingDistance with the following:

SELECT drv.*, rn.geom_way, rn.x1, rn.y1 FROM pgr_drivingdistance('
  SELECT rn.id AS id,
  rn.source::int4 AS source,
  rn.target::int4 AS target,
  rn.cost::float8 as COST
  FROM public.na_2po_4pgr rn
  where ST_Intersects(rn.geom_way, st_geomfromtext(''POLYGON((-79.5670159782642 36.3965968758439,-79.5670159782642 36.6695646217025,-79.2319617622083 36.6695646217025,-79.2319617622083 36.3965968758439,-79.5670159782642 36.3965968758439))'', 4326))',
  14232131,
  25,
  false,
  false) as drv
  left join public.na_2po_4pgr rn on drv.id1 = rn.id;

But I receive no results at all.  At first I thought perhaps the start ID needed to be either the source or the target id, and not the id of the record (though I'd expect it's looking for the id specified by the id column in the inner query).  I tried providing the source (and target) value of the starting record as the input, and interestingly, that DID return records, but in a completely different part of the country, and it included records that aren't even in the sub query, so I have no idea how that would even work.


texts: Convert canvas coordinates to screen position

texts:
I have a point, which has map (canvas) coordinates and it would be great to have his screen position. I need this for setting cursor position in the same place, where point is located.

Do you know any way to convert this coordinates?


texts: Using field created in ModelBuilder in same model?

texts:
I have created a model in ArcGIS ModelBuilder and in a point of the process I insert a python tool (New Field) to create new field in a shapefile.

After that I want to use this new field in other tool in the same model. But I cannot use it because it is not created, it is not selectable in the available fields of the shapefile.

This image shows a part of the model where I create the new field and after I want to use this new field as a input in Pivot Table tool.

The New Field toll python code is:

def getParameterInfo(self):
    """Define parameter definitions"""

# First parameter

    param0 = arcpy.Parameter(
        displayName="Shape entrada",
        name="in_features",
        datatype="GPTableView",
            parameterType="Required",
        direction="Input")

    parameters = [param0]
    return parameters

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True
def updateParameters(self, parameters):return
def updateMessages(self, parameters):return

def execute(self, parameters, messages):
    """The source code of the tool."""

#-----------------------------
camp= "GROUP_"
val= "Pix_Value"

    a= parameters[0].valueAsText

arcpy.AddField_management(a, camp, "TEXT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(a, val, "SHORT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
#-----------------------------

with arcpy.da.UpdateCursor(a, ["CODE2012",camp, val]) as cursor:
    for row in cursor:
        a= int(row[0])
        if a < 12100:
                    row[1] = "Residential"
            row[2]= 5
            elif a==12100 or a==12300 or a==12400:
                    row[1] = "Industrial"
            row[2]= 4
            elif a>12100 and a<12300:
                    row[1] = "Highway"
            row[2]= 6
            elif a == 21000:
                row[1] = "AnnualCrop"
            row[2]= 1
        elif a == 22000 or a== 25000:
                row[1] = "PermanentCrop"
            row[2]=2
        elif a == 23000:
                row[1] = "Pasture"
            row[2]=3
        elif a == 31000:
                row[1] = "Forest"
            row[2]=8
        elif a == 32000:
                row[1] = "HerbaceousVegetation"
            row[2]=9
        elif a == 50000:
                row[1] = "Water"
            row[2]=7
        else:
            row[1]= "No_clasified"
            row[2]=0

        cursor.updateRow(row)
    del a
    del cursor`


texts:
There are two different ways the polygon below can be represented:

In ESRI-land, the polygon would be modeled as a single exterior ring that loops back on itself.
In the OGC world, it would be modeled as a polygon with an exterior ring and an interior ring.

Does the GeoJSON format require that either one of these approaches be used over the other?


texts: ArcGIS 10 exports losing precision with large numbers

texts:
I wonder if anyone else has run into this problem; and I'm baffled as to its reason and solution.

Whenever I do an export of features that have large-number attributes (over 14 digits) to a shapefile, they seem to lose precision. ArcGIS presents them in a default format of e-notation if they are large, but it seems to lose the last few digits of the number, for example:

728400000032249 becomes 728400000032000

Any ideas?


texts:
I'm working on a model that will create different buffers for different selections. So I'm using Make Feature Layer and adding an expression to select a subset of features and then using the buffer tool on that subset. But there may be times when a selection doesn't return any features, in which case I wouldn't want to run the buffer. Is there a way in ModelBuilder to test whether or not any features are selected and only run the subsequent step if there are?

Alternatively I could export to a Python script and handle this testing there but I'm not quite sure how to accomplish this in Python either.


texts: How to add elevation data to tif to create DEM (SRTM satellite data)

texts:
I have a .tif created from the SRTM satellite data. The tif should already be DEM data as it was collected from the SRTM satellite.

When I open the tif in ArcGIS, I don't see any height or z data, nor any way to add it.
What am I missing here? I tried to add XY data as you would with a shapefile, but that option is not available for this tif. I tried "Project Raster" and "Define projection" but again, there is no option to project to a vertical coordinate system.

Where is the Z data in this tif??

The tif was downloaded with kml files--which I think are "hillside shade" --is the Z data somehow in the kml layer?

Edit:
This is a .tif from the SRTM dataset. It's distributed by USGS.
Further information about it:

"SRTM GL1 (30m Ellipsoidal) is a version of the popular Shuttle Radar Topography Mission dataset where elevation values are WGS84 ellipsoidal height as opposed to the standard orthometric, or geoid-referenced, elevation. OpenTopography constructed this dataset by subtracting the EGM96 geoid model from the standard SRTM GL1 data. This dataset was produced primarily to support the InSAR community and their need for ellipsoidal referenced DEMs for terrain correction"

XY Coordinate system: WGS 1984
Vertical Coordinate system:  there is none

Raster information:

Number of bands: 1

Format: TIFF

Statistics

Band 1:
Min -56

Max 22768

Mean 208

Std dev 373


texts: Clip error with QGIS - idles without finishing

texts:
I need help troubleshooting a problem using the Clip tool with QGIS 2.8.1 Wien.

When I try to use the tool, it just sits there idling and never finishes. I've let this run for hours and there hasn't been any change in the progress meter.

I've checked that all of the support libraries (numpy, gdal, python, etc) are up to date - they are - but it doesn't make any difference.

I've made sure that the CRS is the same for both layers.

I've reviewed the Console output on my mac but don't see any errors appearing (maybe there's another log file to check for QGIS?).

Here is a screen shot of what I'm seeing (note that the name for the "Output shapefile" disappears after I click the OK button):


texts: Creating Interactive Map Using Something Other Than Google?

texts:
I have a Roleplay Forum, and I would like to add a map of my Roleplay World to the website.

I would like my map to be clickable, and have each location take you to a place on the forum. I don't know coding outside of basic basic phpBB prosilver.

I want the map to be an image, a large one, where you can scroll around the element like a map and then click on a city, location or land and it takes you there on the forum.


texts: Using if statement with ArcPy and GetCount?

texts:
I am looking to get a count on the polygonNieghbor and if it has polygons in it then I want to move to the next step MakeFeatureLayer. and If it doesn't have any features I would like it to finish.

Is there any other way to do this?

arcpy.GetCount_management(PolygonNeighbor_TableSelect)
arcpy.AddMessage(arcpy.GetMessages())

if arcpy.GetCount_management(PolygonNeighbor_TableSelect) > 0:
    arcpy.MakeFeatureLayer_management(MapUnitPolys, inFeatures_lyr)
 else
     print "done"


texts: Avoiding greyed-out rows in ArcGIS batch processing?

texts:
I'm using the batch processing ability in ArcMap 10.6.1 to run a ModelBuilder tool using Zonal Statistics As Table and Add Fields on a bunch of rasters. When I try to specify where to save the output tables, all but the top line of the batch is greyed out. This is frustrating because when I start filling out all the other fields of the batch, the rows are not greyed out, but by the time I finish they often are. I lose all the time I put into filling out all the fields for the whole batch, only to be able to run the tool for the first raster.

While I am using a tool from ModelBuilder, I have had this happen when batch processing the existing tools too, including Zonal Statistics as Table.

Does anyone know how this happens and is there a way to avoid it?

I do not see a pattern for what causes it to go grey versus keep working for all the rows.


texts: Overflow when using projwin with gdal_translate

texts:
Given the following dataset:

Driver: GTiff/GeoTIFF
Files: eudem_dem_3035_europe.tif
Size is 240000, 200000
Coordinate System is:
PROJCS["ETRS89 / LAEA Europe",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",52],
    PARAMETER["longitude_of_center",10],
    PARAMETER["false_easting",4321000],
    PARAMETER["false_northing",3210000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","3035"]]
Origin = (2000000.000000000000000,6000000.000000000000000)
Pixel Size = (25.000000000000000,-25.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 2000000.000, 6000000.000) ( 51d50'39.04"W, 66d46'12.03"N)
Lower Left  ( 2000000.000, 1000000.000) ( 13d41' 3.88"W, 28d46'45.49"N)
Upper Right ( 8000000.000, 6000000.000) ( 86d 0'31.83"E, 56d26'45.47"N)
Lower Right ( 8000000.000, 1000000.000) ( 46d30'12.69"E, 24d 4'54.00"N)
Center      ( 5000000.000, 3500000.000) ( 20d26' 6.96"E, 54d 9'53.47"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Min=-226.887 Max=5111.089
  Minimum=-226.887, Maximum=5111.089, Mean=341.729, StdDev=493.987
  NoData Value=nan
  Metadata:
    STATISTICS_MAXIMUM=5111.0888671875
    STATISTICS_MEAN=341.72920322335
    STATISTICS_MINIMUM=-226.8865814209
    STATISTICS_STDDEV=493.98749395974

when using gdal_translate to extract something from lower right corner:

-projwin 7999000 1002025 8000000 1002000

an error is thrown

Input file size is 240000, 200000
0ERROR 1: TIFFFillTile:Read error at row 4294967295, col 4294967295; got 0 bytes, expected 1534
ERROR 1: TIFFReadEncodedTile() failed.

ERROR 1: eudem_dem_3035_europe.tif, band 1: IReadBlock failed at X offset 937, Y offset 780
ERROR 1: GetBlockRef failed at X block offset 937, Y block offset 780

It works if extracting something where y>4915900.000, and seeing that it tried to read from 4294967295 which is 2^32 = (4294967295+1)

Anyone with suggestions to what can be done? I am using geodata/gdal for running the commands on docker.

with debug on

GDAL: GDALOpen(eudem_dem_3035_europe.tif, this=0x1364b20) succeeds as GTiff.
Input file size is 240000, 200000
0GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: GDAL_CACHEMAX = 348 MB
GDAL: GDALDatasetCopyWholeRaster(): 330*256 swaths, bInterleave=0
ERROR 1: TIFFFillTile:Read error at row 4294967295, col 4294967295; got 0 bytes, expected 159448
ERROR 1: TIFFReadEncodedTile() failed.

ERROR 1: eudem_dem_3035_europe.tif, band 1: IReadBlock failed at X offset 237, Y offset 417
GDAL: GDALClose(eudem_dem_3035_europe.tif, this=0x1364b20)


texts: How can I maintain MapInfo GRC file attributes in QGIS/GDAL?

texts:
I am trying to manage a classified MapInfo grid file in QGIS - I'm wondering if anyone has any experience with maintaining the attributes of the classified grid?
QGIS seems to assign them just a number (Grid value 51, 52, 53 instead of site1234_sector1, site1234_sector3).

It does seem to recognize that the different attribute values exist, and lists them in the metadata, but I don't see any simple way to pull that table linking the new grid values to the old grid attributes.

Any advice would be appreciated. Thanks!

EDIT: Current system that works involves copying the attribute list out of the metadata tab and pasting it into excel. This works, but clearly isn't ideal.


texts:
I have two feature classes.

Point feature of hand-sampled GPS locations
Line feature of terrace walls digitized from aerial photograph.

Because of the less than perfect ways I generated these features, I need to allocate each GPS point to the nearest terrace wall from which I sampled it.

I defined the Input as the Point Feature, and the near feture as the terrace lines.

But the process fails to generate results.
If it's a Near Table, then it's an empty one, and when I use the Near Tool itself on the point feature, it adds two fields with -1 in them.

When I run the tool on the features in an opposite way, it does generate results, and for every terrace walls gives me the nearest point - obviously, not what I'm looking for.

Any idea why the process fails?


texts: edit openlayers feature attribute and save in geoserver

texts:
I have a layer in geoserver and I want to get this layer with WMS in OpenLayers.

population = new OpenLayers.Layer.WMS("population",
    "http://localhost:8080/geoserver/world/wms",
    {'layers': 'population', transparent: true, format: 'image/png'},
    {isBaseLayer: false}
);

I get the property of feature with OpenLayers.Control.WMSGetFeatureInfo.

featureInfo = new OpenLayers.Control.WMSGetFeatureInfo({
    queryVisible: true,
    layers: [population],
    output: 'features',
    infoFormat : 'application/vnd.ogc.gml',
    format: new OpenLayers.Format.GML,
    eventListeners: {
        getfeatureinfo: showInfo
    }
});

I want to edit the property of some features.

How can i do?


texts: Connect to PostgreSQL from ArcMap 10.2 failed

texts:
I have a database server PostgreSQL and Postgis. I want to connect to database from ArcMap 10.2 , I tried Add Database Connection from arcmap menü.

I typed my connection informations and clicked database selection box, an error occured as following. My connection info is true. I can connect from pgAdmin to database.

I did not install ArcSDE to database server. I read before that I do not need SDE after 10.0 version.


texts: 3D to 2D from point cloud using a georeferenced orthophoto?

texts:
I have 3D point cloud (laser scanning) and an orthophoto, geo-referenced by a TIFF World File (.tfw). I want to project every 3D point of the point cloud into the 2D image.

Here is an example of a .TFW file:

2.4384
0.0000
0.0000
-2.4384
441794.4342
5094101.4520

These are the 6 coefficients of a three-by-three transformation matrix. They indicate that a pixel is 2.4384 meters square, and that the center of the upper left pixel has an easting (x coordinate) of 441794.4342 and a northing of 5094101.4520. See here.

So I can get the X,Y of the point cloud (supposing that the Z is 0):

x_cloud = 2.4384 * column + 0.0 * row + 441794.4342
y_cloud = (-2.4384) * row + 0.0 * column + 5094101.4520

How do I find (X_pixel, Y_pixel) from a 3D point using the TFW file?

Solution

I found a solution for my problem, it was simple. let's see our TIFF World File

2.4384 -> A
0.0000 -> B
0.0000 -> C
-2.4384 -> D
441794.4342 -> E
5094101.4520 ->F

here is my pseudo code:

forEach point cloud  {

    x_cloud = currentX_cloud;
    y_cloud = currentY_cloud;
    z_cloud = currentZ_cloud;

    x_pixel = (x_cloud - E) / A;
    y_pixel = img.rows - (y_cloud - F + A * img.rows) / A;

//accessing to pixel values

    R = img.at<cv::Vec3b>(y_pixel, x_pixel)[0];
    G = img.at<cv::Vec3b>(y_pixel, x_pixel)[1];
    B = img.at<cv::Vec3b>(y_pixel, x_pixel)[2];

}

(I'm using OpenCV to manipulate the orthophoto and LidarFormat to deal with the point cloud data (LAS file))


texts:
I'm confused because of the representation of mbtiles files in QGIS. It seems that the whole grid is shifted by about 3 m to the south.

So far I can observe this behavior only in the context of QGIS and mbtiles.
Maybe somebody comes up with an idea?

For reverse engineering:
Load Geotiff in Tilemill and export as mbtiles. Load both files in QGIS.

On the search for the cause I have

generated the mbtiles file with geoscript-groovy: same result -> Tilemill does not seem to be the cause.
taken an other GeoTIF -> same result: shifted 3m south.
loaded the mbtiles-file in MapProxy, published it as WMS and loaded it in QGIS -> fits perfectly: the mbtiles file itself it seems so not to be the cause.