texts: How to rotate an embedded web map (on pc)?

texts:
I was looking through a couple of maps I made, the maps were created by serving VectorTiles through TileStache and using Leaflet and Tangram (on respective pages) to render them. One of the reasons I wanted to use a Vector Tile map was because Vector Tiles are  rendered on the client end (as opposed to the per-rendered raster tiles we get) so the styling is applied by the client, which means that the labels and icons applied on them will not be static (they should realign when the map orientation is changed). I wanted to test this by changing my map orientation, and realized that I can't rotate my map on the pc. I looked into it and discovered that  Google Maps can't be rotated as well.

Is there any way I can rotate the maps I have embedded? Is there a way to somehow change the orientation of Google Maps (or any embedded web map) on a PC?

Update: So OpenLayers 3 supports rotation but I can't embed Tangram in OL3, is there a way to rotate a Leaflet map?


texts: reduceRegions for multiple images at once?

texts:
Is it possible to get statistics for multiple images at once using reduceRegions() function?

My case: I need to get statistics from multiple images (elev, slope, tree cover loss) for multiple region (grid)

// Add reducer output to the Features in the collection.
var gridmeans = elev.reduceRegions({
  collection: grid,
  reducer: ee.Reducer.mean(),
  scale: 30,
});

print(gridmeans);

Export.table.toDrive({
  collection: gridmeans,
  description:'biggridmeans',
  fileFormat: 'csv'
});


texts: Character encoding support in geodatabases and shapefiles

texts:
I've got several geodatabases that include feature classes with Greek letters in a lot of the attributes.  When I try exporting a feature class as a shapefile from ArcCatalog, the attributes are butchered in the shapefile data, some kind of character encoding problem (they look like this in the shape: ??etr??e?).  The same thing happens when I use ogr2ogr in FWtools to convert layers from the MDB to KML, shp, etc.

Does anyone have experience trying to deal with encoding formats across GIS data formats?

The real goal here is to get some data from these Esri geodatabases into a Postgres/PostGIS database, but broken encoding isn't going to work.  I was going to export from the geoDBs to shapefiles, then load them in with shp2pgsql.  Is that the easiest path to get there?


texts: Connecting to Oracle database with ArcGIS for Desktop?

texts:
I installed ArcGIS 10.2.2 ,oracle 11g R2 personal edition and it's client.I try to connect to the database using database connection in ArcMap but the error

  Failed to connect to the specified server. Do you want to continue?
  Cannot connect to database because the database client software failed
  to load. Be sure the database client software is installed and
  configure correctly.

raised when i want to connect . How to solve the problem ?


texts: What is the future of CartoCSS (i.e. does it have one)?

texts:
MapBox recently announced that they would be moving away from CartoCSS towards a completely visual interface.  Apparently the new style sheets that underlay that UI will be JSON based, so despite the fact that the new interface is a point and click system one of the developers that worked on this project believes that it will ultimately be better for coders than the previous system:

  the
  new editor https://t.co/i7uy2V8YGL is even
  better for coders: the GL styles can be modified as JSON in editors,
  and dynamically in JS— Tom MacWright (@tmcw) July 29,
  2015

It's worth noting however (and hopefully I have this right) that these JSON stylesheets are rendered by WebGL, while CartoCSS is used to create XML that is rendered by Mapnik.  So while one is a replacement for the other, at least in some capacity, they don't do exactly the same thing.

So this leads me to wonder where does that leave CartoCSS?

Many of the predominant map tile styles use the language at this time, perhaps most notably OpenStreetMap's default style.  CartoDB uses it as well and as one of their developers noted in the tweet below they plan to stick with it for the time being, but it also seems to imply that they will move away from it at some point.

  @Anonymaps @tmcw we however are sticking with
  it a lilttle longer!— Francisco Dans (@fdans) July 29,
  2015

On a more recent front there's a project that called kosmtik which I found through this thread.  This tool reads the same .mss (cartocss) and .mml (json) files as TileMill and renders with Mapnik so it does provide an alternative to TileMill which is no longer being actively developed.

The ultimate source of this question is that I undertook a pretty complex project a month or two ago in TileMill, and am perhaps 25% to 50% towards completing it.  But with this announcement among other things I'm beginning to wonder if I'm hitching myself to a medium that won't be maintained/developed going forward.

So my question is, what other evidence exists that speaks to level of use of CartoCSS (or lack thereof) going forward by the GIS community at large?


texts:
I've this problem when I try to use an URL of a map WMS for example

Prueba de conexión fallida [Test connection failed]: Server returned HTTP response code: 401 for URL: http://focosdecalor.conae.gov.ar/geoserver/FocosDeCalor/wms?REQUEST=GetCapabilities&VERSION=1.1.0&SERVICE=WMS

what is this?


texts: Obtaining date data from layer name in QGIS

texts:
I have a column with name of the satellites image background (Landsat, Sentinel), in this name I have a date information into the text.
How do I extract the date data from this text?


texts:
We are in the process of getting our property boundaries corrected and a side effect of this is that we need to move any assets that have been drawn relative to the property boundaries, so they are still inline correctly.
The company that does the property corrections gives us shift vectors and we normally run all our layers though FME using the vectors to correct the locations.
I was wondering if there is anything like this in GRASS?


texts: FME Server API - Retrieve Output file as a URL

texts: I have a FME Server web workspace converting a file from SHP to KML, once it runs it outputs a kml file, I need to be able to take this file and drop into my web API google Map. Is there a way to retrieve this temporary file and display in my web map?


texts:
I have a geodatabase (ESRI/PostgreSQL) I'm working in an application for digitizing :
in my country we have 4 coordinate system(projection), so I'm looking for a way to switch projection or any other idea or set projection parameter to all my feature classes.


texts:
This is what my screen looks like now.

The yellow outline of the US is a shapefile and so are the blue dots. I am using one of the default basemaps in the background. I want center the map so that the entire US is visible. I've tried changing the scale but it hasn't fixed the issue. I can't simply drag the map because then it will be off the page I want to print.


texts: MapGuide and Google Maps Coordinate System problems

texts:
I'm trying to do a simple MapGuide OS + Fusion overlaying Google Maps and I can't figure out why the coordinate system are off. I downloaded populated places from natural earth.

After Googling, it says that Google Maps uses WGS84 and natural earth is WGS84. Here's what I did:

I create a map definition from the natural earth data. Then I created a Flexible Layout with that map definition to view it with fusion. I checked the coordinate system of the map. It says "GCS_WGS_1984"

When I display the map, it is incorrect. (The points are way to small and do not correctly align)

How do I go about trouble shooting this? I'm brand new to GIS and I think it has to do with the coordinate systems, but as far as I can tell they should match.


texts: Converting polygon with latitude and longitude to x y plane

texts:
I have a polygon (shape of a building) containing coordinates as latitude and longitude. I need to convert them into an x-y plane. I have used various methods such as:
const R = 6.371;
const x = R * Math.cos(lat) * Math.cos(lon);
const y = R * Math.cos(lat) * Math.sin(lon);

However, all methods do not scale up the coordinates in the x-y plane. The result ends up being a tremendously tiny polygon.


texts: Are Database Connections the best way to access Feature Classes?

texts:
I wrote some code on my machine, debugged, tested... it works great. Now, a coworker wants to use my scripts. The problem, I used database connections to reference the database and feature classes and everyone names their database connections in a slightly different way. Is there a better way to reference the SDE database and features classes?

Here's how I reference the objects.

outputFC = "Database Connections/XXX_YYY_PROD.sde/Schema.FeatureClass"

I did declare variables to contain the connection string but nonetheless someone still needs to change the script, even if it's only in one place. Is there a way to reference the database directly?


texts:
I have a point shapefile on which I would like to use the "Add geometry attributes: tool in order to add x and y coordinates using a GCS that I am interested it. When I run this tool I get the runtime error:

  "The field is not nullable (Point_M)".

I don't understand why I am getting this error. AddGeometryAttributes is supposed to add this field on its own, I have no role in the field Point_M. In fact it doesn't even exist before I run the tool. This is one of the outputs of AddGeometryAttributes.

Has anyone had this problem before?


texts: How to fix the "The inserted Ring is not contained in a feature"?

texts:
I am trying to add a ring in a feature in a shapefile. This feature is a single part feature but has already a few rings. When I try to add the new ring, I have the error message : "The inserted Ring is not contained in a feature". I tried to modify my snapping options thinking that the problem came from this setting, but I didn't have any luck.

Does anybody know how I can fix this problem?

I use QGis 1.8.0-Lisboa and Windows 7.

Thanks


texts: Use mapview r package to display vector data with outline only

texts:
I'm creating a map in r using the mapview package. I have a raster layer with a vector layer displayed over it. I want to view the raster data completely. How do you set vector data to be completely transparent (i.e. just an outline)?

A reproducible example from their vignette:

m1 <- mapView(franconia, col.regions = "red")
m2 <- mapView(breweries)
### add two mapview objects
m1 + m2
'+'(m2, m1)
### add layers to a mapview object
if (interactive()) {
library(plainview)
m1 + breweries + plainview::poppendorf[[4]]
}

In this example, I would want breweries to be completely transparent. However changing col.regions = 'transparent'  for breweries isn't what I'm looking for. It's still opaque as seen when clicking on the zoom button (plainview::poppendorf[[4]]) on the lower right hand side.


texts: Using ArcMap I want to add new Fields into a gdb using Python

texts:
I'm new learning as I go, I have read some of the Answers posted, and I've gotten this far, I was able to add a new field to a GBD, I wanted to have the Length be 50, but it created a field with 255, so I didn't get something quite right  the code I used  below  the feature class is wax_ru the Field name is HN_AW  Text is the Type of Field, and 50 is the Default Value I wanted, but as I stated it created the field with 255

import arcpy
from arcpy import env
env.workspace = "C:\Users\angelaw\OneDrive - TomTom\Resources\Australia\PSMA_MAY_2018\PSMA_Deltas_test\PSMA_MAY_deltas.gdb"

arcpy.AddField_management ("wax_ru", "HN_AW","TEXT",50)


texts: PostgreSQL: Spatial join points to polygon, listing point attributes where multiple are joined

texts:
I have a large polygon table and a large point table in a postgis enabled PostgreSQL database.
Each point has a unique reference number (double type). I want to spatial join the points to the polygons, but where multiple points are found within a polygon, I want ALL the reference numbers to be carried across into a new column in the polygon table (a string column called ref_numb_list). Multiple reference numbers would be delimited with commas.

This functionality is available in ArcGIS spatial joins where you're able to set a 'join' merge rule for an output column. But I need to do this in PostgreSQL, via pgAdmin.


texts: Merge String fields in ModelBuilder by undefined ArcGIS user?

texts:
I'm using ArcGIS 10. I have a doubt, because I already have a merged shapefile, with two string Fields, some Null and some filled (comes from the merge tool).

So, now I want to get the values from the two string fields into a NEW string Field but with ModelBuilder.

So, I do know how to do it manually and get these values with the FIELD CALCULATOR, because I could choose the field names, but how could I get this prepared for an user with ModelBuilder?

VB Code:

Dim X
if IsNull( [ETIQUETA] ) then
   X = [Data_Hamlet_Galicia_POI_NAME]
else
   X = [ETIQUETA]
end if
__esri_field_calculator_splitter__
X

The FIELDs above may not have the same name, so it will not be right to let this code into the "Field Calculator" Tool within ModelBuilder.

Any suggestions? Here is my model.

--------   E   D   I    T    ------

Well, @Hornbydd , I've tried but fails, maybe I'm doing something wrong. Here the new model and the error log:


texts: Need advice for localized survey base point

texts:
I have a set of 4 control points used for H and V localization. There is a base point located far to the SW of the property. I did not create this myself, but took it off of an existing Pocket 3-D project from a past employee. I didn't know any other way to find a base point b/c the company is too far out in the country for existing survey marks. The data from points and lines I create in this localization all end up inside the control points where they should be as I set up the base station over the control points 1-4, and the base control point only has N 5000, E 5000, and Z 1000, so no way to figure out where it really is. However, when I got on the computer version on the program, and looked at the list of control points, it shows the WGS84 locations for the control points. The base's WGS84 location in actually inside the control points, so now I'm very confused. I can't figure out how to get WGS84 coordinates for the points, polylines, ect, and have been trying to shoot it though the program's coordinate calculator, getting results saying I'm in eastern Colorado, ect. I'm in North central Kansas, so that's not even close.

I understand that I may need to just start over, and if I do that, I'll need some serious survey advice. I'm an ArcGIS user, stuck doing my first real taste of surveying in the context of construction grade, sub-centimeter level, products (Topcon, Pocket 3-D, 3D Office). I have no training in this really, and there's no one at this company for me to talk to, as I'm basically my own department. Please help! I got in a bit deeper than I expected!

BASE is N 5000', E 5000', E 1000'
    * N38.561866674', W99.212131227', H 2081.0391'
CP 1 is N 7443.5450', E 6320.5000', E 2065.9120'
    * N38.561000637', W99.214876310', H 2089.3173'
CP 2 is N 10039.8030', E 6336.8390', E 2064.4450'
    * N38.563566563', W99.214856097', H 2087.8405'
CP 3 is N 10040.2260', E 8942.3830', E 2080.0470'
    * N38.563567238', W99.211558647', H 2103.4836'
CP 4 is N 7445.9530', E 8922.2990', E 2032.7590'
    * N38.561003304', W99.211584006', E 2056.1248'

the other points in the screen shot were used to get closer to the area I was surveying with the base receiver, but were not used for Horizontal or Vertical localization. Also, the error values on the control points are now higher than they were when I started, is that normal?

This is what's coming up on the program. It sits where it should on Google Earth, maybe slightly skewed, but the lat and long is different from what is in your picture from QGIS.


texts: How to populate fields for line segment with value from polygon to left and right of it?

texts:
I have a countywide roads dataset that is used by an emergency dispatch system.  There are two fields, "Comm_RT" and "Comm_LT", that need to be populated by the name of the community polygon they reside in.  I am having trouble finding a way to populate these two fields when the road borders two different community polygons and thus requires that the "Comm_RT" and "Comm_LT" fields be different.  I have started by making selections in ArcMap 10.2 of all the roads that intersect two different polygons. Here is an image of a road segment that would have two different values for "Comm_RT" and "Comm_LT":
Differing Comm Codes http://slu.s3.amazonaws.com/comcode.PNG

Is there a way to perform a select by location query or a tool in the Network Analyst extension that could be used to solve this?


texts: ArcGIS Extract by Attributes Error: "A column was specific that does not exist."

texts:
I have a single band DEM raster. I am attempting to use the tool "Extract by Attributes" but am getting the error:

"A column was specified that does not exist"

My raster is a single band raster. My understanding is that thus the only column name is "VALUE". I am trying to extract a raster where the value is between a certain range. To ensure I have entered the SQL query correctly I have used the SQL Query creator which is built into the "Extract by Attributes" tool. Here is an example of my inputs with the query shown:

[

What is the error in my approach?


texts: Creating spatial join between points and polygons in GeoPandas?

texts:
I want to create a simple spatial join between points and polygons using GeoPandas but I think GeoPandas has bug?

GeoPandas code:

from geopandas import gpd
import geopandas
points = geopandas.GeoDataFrame.from_file('points.shp') # or geojson etc
polys = geopandas.GeoDataFrame.from_file('polygons.shp')
pointInPoly = gpd.sjoin(points, polys, how='left',op='within')

Error :

Traceback (most recent call last):
  File "/home/username/testshapely/sumpointsinsidepolygon/testgeo.py", line 7, in <module>
    pointInPoly = gpd.sjoin(points, polys, how='left',op='within')
  File "/usr/local/lib/python2.7/dist-packages/geopandas/tools/sjoin.py", line 57, in sjoin
    r_idx = np.concatenate(idxmatch.values)
ValueError: need at least one array to concatenate

And if I change the imports with the some code:

import geopandas
import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame, read_file
from geopandas.tools import sjoin
from shapely.geometry import Point, mapping,shape
import pandas as gpd

I get the following error:

    pointInPoly = gpd.sjoin(points, polys, how='left',op='within')
AttributeError: 'module' object has no attribute 'sjoin'

Any ideas why?


texts: Sorting cursor made of Network Analyst outputs?

texts:
I have solved OD Cost Matrix analysis and I'm trying to access the results using da.SearchCursor with SQL sorting. The problem is that the result is not sorted at all.

lines = arcpy.da.SearchCursor('OD_potential\Lines',['DestinationID','Total_'+impedance],None,None,False,(None,'ORDER BY DestinationID ASC'))

This approach works with every other layer.

Previously I used old Cursor and it worked fine, but it was much slower.

lines = arcpy.SearchCursor('OD_potential/Lines','','','DestinationID; Total_'+impedance,'DestinationID A')

I know there is a workaround using arcpy.Sort_management before the Cursor but then the runtime is same as using the old Cursor.

Any suggestions?


texts: How can all polylines within a polygon be selected and then operated on to determine the average value from an attribute field?

texts:
I have two files illustrated wonderfully below: 1) Polyline (shown in red), 2) Polygon (shown in black).

The polyline file contains a numeric attribute (elevation in meters).

I would like to select all polylines that are contained within (and intersect) each polygon grid cell. I would then like to count the number of lines inside each cell, sum their elevations (from the attribute table) and report the average elevation of polylines in each cell. The lines can belong to more than one cell.

Aside from doing this manually using selections, which is not possible given the number of cells/lines, how can I accomplish this task?

In the end I would simply like a list of the average elevations for each grid cell - ideally this could all be done in the GIS, but exporting to a spreadsheet and performing the final calculations (average) is workable.


texts:
Has anyone been able to keep multiple popups open at the same time -- Leaflet output -- from the latest version of QGIS2Web?

I have  a leaflet map generated from QGIS2Web that has several different point layers and keeping multiple popups would be highly beneficial.

In the pre-Leaflet 1.0 version of the QGIS2Web plugin I was able to keep multiple popups open by inserting this code in the HTML sheet:

L.Map = L.Map.extend({
openPopup: function(popup) {
            this._popup = popup;

    return this.addLayer(popup).fire('popupopen', {
        popup: this._popup
    });
}

});

I tried slotting this into the leaflet 1.0 version output and I can't say that it worked. I've been poking through the HTML file to try and find a solution but I'm coming up short (can't say I'm the greatest with HTML/JavaScript). Has anyone had success with this?


texts: Change projections depending on the zoom level

texts:
Is it possible in ArcGis (from a technology/performance/… point of view) to switch projections depending on the zoom level?

For example, Is it possible when zooming at some extent to the poles, my projection to be switched from the "Miller projection"  to "stereographic projection". For polar areas, a stereographic projection is best.


texts: PyQGIS QgsProcessingAlgorithm: no output layer when cancelled

texts:
For one plugin I maintain, I can't figure out how to return output layer when the user hits 'Cancel'. Basically I'm looping over a set of results in processAlgorithm which populates the output layer and I'd expect to have the output layer returned with all results calculated until 'Cancel' was clicked.

A minimal example:

def processAlgorithm(self, parameters, context, feedback):
    (sink, self.dest_id) = self.parameterAsSink(parameters, 'Output', context,
                                                list_of_fields,
                                                QgsWkbTypes.Polygon,
                                                CRS)

    for result in results:
        if feedback.isCanceled():
            break

        # Do smth smart

        sink.addFeature(result)

    return {'Output': self.dest_id}

However, no layer is added to the map when 'Cancel' is hit mid-way through processing.

EDIT
Just to clarify, I'd have expected that if feedback.isCanceled(): break would yield return {'Output': self.dest_id}, i.e. the output layer with the added features until feedback.isCanceled() is called.


texts: Is there no EPSG code to use for Behrmann equal-area projection in ogr2ogr?

texts:
I'm trying to use ogr2ogr to reproject some shapefiles. The projection I wish to use is Behrmann equal-area. I cannot find an EPSG code to use for the '-t_srs' part of the syntax. Is there any way around this? The best I can find is this - http://spatialreference.org/ref/esri/54017/

Any help much appreciated!

Update:

for file in $INPUT/*.shp ; do
ogr2ogr -t_srs PROJCS[\"World_Behrmann\",GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984  \",SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Behrmann\"],PARAMETER[\"False_Easting\",0],PARAMETER[\"False_Northing\",0],PARAMETER[\"Central_Meridian\",0],UNIT[\"Meter\",1],AUTHORITY[\"EPSG\",\"54017\"]] $OUTPUT/$file $file
done

This is now what I have. Apart from inverted commas, i'm not sure what else needs to be escaped in this statement. Any ideas?


texts: Transform data type EWKT and GeoJSON to SRID 900913


texts: Finding applied transformation method and accuracy information when reprojecting layer in QGIS

texts:
When I load two layers with different CRS I'm prompted to choose a transformation method and there is a column indicating the accuracy (see image 1). I would like to know where the information as part of the on-the-fly (OTF) transformation is coming from.
As far as I can tell, there is no such information displayed when I try to reproject a layer into a new CRS which is based on a different ellipsoid, using the 'Reproject Layer' tool. All I seem to be able to do is to choose the target CRS but there is no option to choose a particular method, neither there is the info about which method is applied by QGIS, and what the accuracy implications are (see image 2).
How can I address that particular component shown as part of the OTF with Python?
My research indicates I have to make use of the proj4 library but I don't know how to do it yet.
Is the transformation code (tfm code - see image 1) also part of the library and how to address it?
Background of my question: As part of my coursework I intend to write a plugin which lists possible transformation methods and the accuracy when using the Reproject tool, similar to the OTF transformation. If this isn't possible (I can imagine there is a reason why it is not already there), my alternative goal is to display the information about which method is applied and the accuracy at least.
I have carefully read the related posts on this plattform, however I get the impression they don't help with what I'm after.
I am using QGIS version 3.10.5.


texts: Applying DEM elevation adjustment parameters in QGIS?

texts:
I have a nice DEM, created from UAV data, which I want to process further in QGIS, my platform of first choice. In order to represent the evelation of each pixel correctly, I have to adjust the 'elevation adjustment parameters' like Z-scale and Z-offset. In Globalmapper, I was able to do it. Also in the manual of ArcGIS, it is discussed (http://resources.arcgis.com/en/help/main/10.1/index.html#//009t0000003t000000).

I do, however, not encounter the way how to do it in QGIS.

Does have anyone a clue or advice on how to do that in QGIS?


texts: QGIS Google satellite/aerial suddenly not displaying

texts:
QGIS 3.10 on Linux Mint 20.
Just got back to work after the Christmas break and a QGIS project which was displaying Google aerials through XYZ tiles perfectly beforehand, is now blank.  (as is Google Maps but G Roads, Hybrid & Terrain are working as are Bing & OSM. )
I'm assuming Google have changed something in the URL, but all the references I've found are old and don't work.
Any insight into this or where I could find the latest URL?


texts:
I have a single band raster layer, namely the USGS Global Mountain Explorer Home K3 dataset.

It is a 8 bit unsigned integer GeoTIFF of 200 Mb, with 4 grey values each corresponding to a mountain height.

What is the correct way to split this layer into 4 different layers so that I can vectorize them separately?

If I mess around with transparency and then I try to save the rendered images the results are huge, even with high compression.


texts: Exporting selected feature to Google Earth in QGIS?

texts: Does anyone know if there is a plugin for QGIS that will export features to Google Earth similar to the one for MapInfo?


texts: Legend only to display what is shown in Print Composer QGIS 1.8.4

texts:
I am using the 625k bedrock data from BGS. I have a range of site across scotland.  I have then clipped the file down to within a 70Km radius of each site, a lot of over lapping areas with disintegrated boundaries ends up with most of scotland covered.

However, i have categorised the style of the bedrock data by name.

When i place the legend in there seems to be an awful amount of classes.

Bar going through each layer and ascertaining if it is displayed or not is there a way to display only the parts of the layer you can see in the legend.

Cannot show images due to client confidentiality.


texts: Proper styling of help-text in custom functions (@qgsfunction) - html-code

texts:
I would like to have a QGIS-like-styling of the help-text of custom functions.

I.e. the html-code
text
 creates the green headline.
Are there theme fonts or other definitions (colors...) for the other elements of the help-text known (marked by <?>...</?> in the picture below)?

Alternatively - where can I find the files with the core-expressions to look there for the definitions of help-text in my qgis-installation-directory?


texts: Mexican "Área Geoestadistica Básica" OSM admin level

texts:
I found shapefiles that split urban areas in Mexico into distinctions called Área Geoestadistica Básica (AGEB) from Mexico's geographic information site, INEGI

This is what the divisions of Mexico city look like:

The link to download the shapefiles is here (32 mb)

I know the following mapping of geographic distinctions exists within Open StreetMap in Mexico:
 (from wiki.openstreetmap.org/wiki/WikiProject_Mexico)

I was wondering if anyone knows if there is any OSM mapping for these Basic Geostatistical areas?
If not, are they just like census tracts/blocks in the US that serve purposes for information gathering, but are not used for administrative purposes?

My apologies if this is more of a basic translation questions - I've looked around for anything in English describing what they are and was wondering if anyone had any more specific information about the nature of these divisions?


texts: Reshaping object in arcpy.da.UpdateCursor

texts:
I want to make a kind of Resape script tool for highlited polygons.
For example there are two polygons from different *.gdb files that share a same boundary and they need to be merged.
I read coordinates form the first shape then from the second and made a list of coordinates. Then united them into one list made an array with arcpy.Point(x,y) features. Finally, did arcpy.da.InsertCursor thing and got a new polygon with geometry of two polygons.
And is there any other way to do it with arcpy.da.UpdateCursor?
For instance, using something like

with arcpy.da.UpdateCursor(polygon_fc, "SHAPE@XY") as uc:
    for row in uc:
        p_obj = uc.getPart(0)
        #some stuff
        uc.updateRow(row) - #calculating "SHAPE" field to bring new coordinates from list of new points

I know that there is an easy way to make it manually through Merge tool, but anyway is it possible to do with a script?


texts: CartoDB JSONP callback in leaflet-search plugin

texts:
long time reader, first time question-asker etc...

So I'm trying to implement a search function on an app I'm building using data from cartodb in leaflet and this leaflet-search plugin: http://labs.easyblog.it/maps/leaflet-search/

As a bit of background, I've got a standard search working fine (calling my geojson using the leaflet-ajax plugin), but the problem is this only includes one field and I need it to pass multiple fields through the search. I know it's possible to do, and I think I've found a method to do it, but something's going wrong and I don't know why.

This example: http://labs.easyblog.it/maps/leaflet-search/examples/jsonp-filtered.html uses a jsonp call from mapquest, then splits the fields (don't ask me how this is useful, but it's just an example), so I'm pretty sure all it needs is to concatenate fields in the "key" var to make them search as one. code from example below

var jsonpurl = 'http://open.mapquestapi.com/nominatim/v1/search.php?q={s}'+
               '&format=json&osm_type=N&limit=100&addressdetails=0',
    jsonpName = 'json_callback';
//third party jsonp service

function filterJSONCall(rawjson) {  //callback that remap fields name
    var json = {},
        key, loc, disp = [];

    for(var i in rawjson)
    {
        disp = rawjson[i].display_name.split(',');

        key = disp[0] +', '+ disp[1];

        loc = L.latLng( rawjson[i].lat, rawjson[i].lon );

        json[ key ]= loc;   //key,value format
    }

    return json;
}

var searchOpts = {
        url: jsonpurl,
        jsonpParam: jsonpName,
        filterJSON: filterJSONCall,
        animateLocation: false,
        markerLocation: true,
        zoom: 10,
        minLength: 2,
        autoType: false
    };

map.addControl( new L.Control.Search(searchOpts) );

In my project I've got the concatenation working fine, but I can't for the life of me get my data into the filter. My code is identical save for the jsonp url and the callback

        var jsonpurl = 'http://stratiform.cartodb.com/api/v2/sql/?q=SELECT%20name,%20ST_AsGeoJSON%28the_geom%29%20FROM%20schools_points%20WHERE%20name%20ILIKE%20%27%25clunes%25%27&format=json',
        jsonpName = 'callback';

I've tried setting '&callback=whatever' on the url, tried using the ST_AsGeoJSON sql function, and I'm almost out of ideas!

Is there a problem with the cartodb jsonp url I'm using? Does it handle callbacks in a different way or something? Do I need to modify it to somehow to match this function? Is it a cross-origin thing?? I guess what I'm asking is why are the mapquest jsonp and cartodb jsonp links working and not working respectively?

I've got work-arounds for this problem, but it's obviously just cleaner to set this stuff up in the app itself. Can anyone help?


texts: How to add links to PDFs in a Shapefile?

texts:
I want an attribute in a shapefile named "Link" to add.

The works are all fine but then I would like to add a link in my directory.

This link accesses a PDF. Now, this PDF is not shown.

Is there any other solution or there is nothing to be done?


texts: Hosting ArcGIS Web Application on AWS S3 as Static Website

texts:
I am trying to host an ArcGIS web application that I built using the Web AppBuilder. I downloaded the zip file of the completed app from Web AppBuilder and uploaded the unzipped contents, without changing anything, to a S3 bucket. I enabled static website hosting and desiginated the index.html and error.html file. I also created a bucket policy to make all of the files public read.
When I go to the endpoint, the blue ArcGIS loading screen appears, but eventually it times out. See below:

I have a feeling it has something to do with the config file, but I don't know what. I played with it for a while, but couldn't figure it out. Does anyone have any ideas?


texts: ArcMap Model Builder Iterator to create selections

texts:
This is a sequel to my question about aggregating polygons in order to create neighborhoods.

I want to automate the geoprocessing of the above task.  When I use Model Builder's "Iterate Feature Selection" iterator I get the below error:

I originally tried to link the "Iterate Feature Selection" iterator directly to the Aggregate Polygon tool.  But as per this thread I inserted the Make Feature layer tool in between.  (Now removed, see comments below)

My problem is I'm unsure how to link the iterator's selection to the input of the Aggregate Polygon tool.  Really, why does the GROUP BY function of "iterate feature selection" iterator not create a complete SQL condition statement??

The following threads have proved un-useful as well:

Iterating in ModelBuilder?

Model builder: select feature based on attribute, iterator?

thanks in advance


texts: On the fly transformation From WGS84 to NAD27 with proj4 code

texts:
I want to use this proj4 code:
+proj=utm +zone=19 +ellps=WGS84 +toclrk66=-3,142,183,0,0,0,0 +units=m +no_defs
or this one
+proj=utm +zone=19 +ellps=WGS84 +tonad27=-3,142,183,0,0,0,0 +units=m +no_defs

...to make on the fly transformations from WGS84 to NAD27, using the geocentric tanslation parameters recommended for Dominican Republic. I typed this codes using the "New CRS" capacity in QGIS 2.0 (Configuration>Custom CRS), but it doesn't accept the transformation parameters. After typing the codes, and pressing OK, the CRS created looks like this:
+proj=utm +zone=19 +ellps=WGS84 +units=m +no_defs

QGIS deletes the geocentric traslation parameters. Maybe I'm violating proj4 codes rules, and that's why I'm asking for help.

I looked up for a solution amongst similar answers on the forum, but couldn't find anything. Sorry if I'm duplicating a question.

Any help would be appreciated.

Jose


texts:
I am trying to translate a .tif raster to an xyz one in Python

The code is this:

import sys
from osgeo import gdal

gdal.AllRegister()
driver = gdal.GetDriverByName('Local Disk')

dataset = gdal.Open('C:\Delft3D\Teste\GDAL\clc_raster.tif')

if dataset is None:
    print "Could not!"
    sys.exit(1)

gdal_translate -of XYZ dataset.tif C:\Delft3D\Teste\GDAL\clc_raster.tif.xyz

dataset = None

But I get the following error.

   gdal_translate -of XYZ dataset.tif C:\Delft3D\Teste\GDAL\clc_raster2.xyz
                        ^
   SyntaxError: invalid syntax

I don't know why this error appears. I used the tips from convert geotiff to simple xyz elevation file?


texts: Deleting rows using UpdateCursor and date field

texts:
I have a large script, in which a part of it wants to check whether a shapefile has attributes in the "Valid_From" date field. If they match the date variables already created, these need to be deleted. Im testing this using a shapefile that already meets the criteria, however the rows are not being deleted.

from dateutil.relativedelta import relativedelta
today = datetime.date.today()
date_after_month = datetime.datetime.today() + relativedelta(months=1)
valid_from = date_after_month.strftime('01/%m/%Y')
valid_to = str(calendar.monthrange(date_after_month.year, date_after_month.month)[1]) + date_after_month.strftime("/%m/%Y")

Now, in the earlier parts of a script, I update a different shapefile using the date variables set above successfully with the following code:

with arcpy.da.UpdateCursor(rtc_shp, ["valid_from", "valid_to", "type"]) as cursor:
    for row in cursor:
        row[0] = valid_from
        row[1] = valid_to
        row[2] = "Offshore"
        cursor.updateRow(row)

However I now have to check whether a seperate shapefile already has this data, and if so, delete the rows. I use the following code, but it does not give the same results. It wont delete the rows. Which is strange, as the same variables are used to update the fields in the above example.

with arcpy.da.UpdateCursor(corporate_rtc_dataset, ["valid_from"]) as cursor:
    for row in cursor:
        if row[0] == valid_from:
            cursor.deleteRow()

print row
2016-10-01 00:00:00

print valid_from
01/10/2016

UPDATE:

I have now tried to change the valid_from variable to match the row, however again, its not deleting the row:

valid_from = date_after_month.strftime('%Y-%m-01 00:00:00')

print valid_from
2016-10-01 00:00:00


texts: Connect ArcMap (Interoperability Tool) to GeoServer

texts:
I have GeoServer (2.8.5) and multiple layers on it. If I want to edit layers in QGIS everything works fine (auth / show layer catalog / load layer / edit layer / save changes).

If I attempt doing the same in ArcGIS (10.2) with the Data Interoperability tool I struggle on getting the capabilities document. The URL is the same but I get following error:

(old image, url's not correct in the image, WFS instead of WMS)

I tried different styles in the URL like /geoserver/wfs? , /geoserver/wfs?request=getcapabilities but nothing seems to work. I've also tried connecting with /geoserver/ows? like mentioned in What's the correct link for adding GeoServer WFS to ArcMap? but that fails too.

Anybody has an idea what I am doing wrong?


texts: Can I make a clickabe image map in geoserver?

texts:
I mean when a region of a map (for example, a country) is clicked to redirect me to another URL, not only to give me some information like fid,name, etc.
You know when a region of a shapefile is clicked below the map is retrieved some information and my question is : would it be possible when the same region is clicked the user to be redirected to another page ?
I am talking about SHAPEFILES not IMAGES(PNG,JPEG,etc) !!
Can this be done?


texts: PostGIS : ST_Contains doesn't work after ST_FlipCoordinates


texts: Tiling triangulated irregular network (TIN)

texts: I want to generate vector tiles for TIN which I created from elevation data. I want to use Tile Map Server (TMS) specification for vector tiles. But I did not understand how tile bounds can be calculated for vector data such as TIN and how quadtree structure can be created for TIN. Any help?


texts: QGIS options menu causes freeze on Ubuntu 16.04, solutions?

texts:
First time asking, I've checked the internet and no luck.  I am unable to edit the program options in QGIS 2.14.3 running on Ubuntu 16.04 (settings -> options).  I've uninstalled all QGIS packages (via purge) and reinstalled, yet the issue still persists.  Regular QGIS operation works, although slow, hence why I need to edit default settings (multipoint shapefile of only a couple megabytes won't render properly).  Does anyone have suggestions for how to do a proper full uninstall and reinstall, or ideas for underlying non-QGIS libraries that could be causing the issue?  Does not appear to be a documented issue online.  The issue persists even when no project is loaded.  QGIS CPU usage goes to 100% and the program has to be killed from command line.

This is the output on the command line both with and without sudo :

~$ sudo qgis
Warning: loading of qgis translation failed [/usr/share/qgis/i18n//qgis_en_US]
Warning: loading of qt translation failed [/usr/share/qt4/translations/qt_en_US]
Warning: QCss::Parser - Failed to load file  "/style.qss"
^C^C^\

~$ qgis
Warning: loading of qgis translation failed [/usr/share/qgis/i18n//qgis_en_US]
Warning: loading of qt translation failed [/usr/share/qt4/translations/qt_en_US]
Warning: QCss::Parser - Failed to load file  "/style.qss"
Warning: QVariantMap DBusMenuExporterDBus::getProperties(int, const QStringList&) const: Condition failed: action
Warning: QVariantMap DBusMenuExporterDBus::getProperties(int, const QStringList&) const: Condition failed: action
Warning: void DBusMenuExporterPrivate::addAction(QAction*, int): Already tracking action "Terms of Service / About" under id 259
Warning: void DBusMenuExporterPrivate::addAction(QAction*, int): Already tracking action "&Toolbox" under id 366
Warning: void DBusMenuExporterPrivate::fillLayoutItem(DBusMenuLayoutItem*, QMenu*, int, int, const QStringList&): No id for action
Warning: void DBusMenuExporterPrivate::fillLayoutItem(DBusMenuLayoutItem*, QMenu*, int, int, const QStringList&): No id for action
^C^\Quit (core dumped)


texts: How to force line dropping in tippecanoe

texts:
I am using tippecanoe to generate tiles for Mapbox.  I am using a geojson dataset composed of lots of lines - I want lines shown at a low zoom to besimplified so it doesnt look too cluttered, but I can't find the right command.  Here are the options I'm using and what I think they are doing:

--drop-rate=2.5 (I think this is the default)
--drop-lines (make dropping of points applicable to lines)
--base-zoom=13 (drop features below zoom 13)
-S 10 (simplify lines)
-z 16 (max zoom=16)
-Z 8 (min zoom 8)

This doesn't seem to work, I get the same number of lines at all zoom levels.  I'd like to know what parameters I should be using to be able to drop lines.
Thanks.


texts:
I am writing a lot of commands in the GRASS GIS console and I get suggestions for auto-complete, i.e. for mapset and layer names.

Unfortunately, this works for layers and mapsets that are already in the location when I first start GRASS GIS. If I create a new layer or mapset it will never show in the auotcomplete suggestions until I restart GRASS GIS.

Is there a way to "refresh" the autocomplete list?


texts: interlis/ili2gpkg: export to XTF with multiple geometries

texts:
I am trying to export a XTF-File from a geopackage with ili2gpkg. The Interlis-Model is KbS_Basis_V1_4/KbS_LV95_V1_4. The data consists of point-features (table "belastete_standorte") from which some also have an additional polygon-geometry (table "polygonstructure"). The data export via

java -jar ili2gpkg.jar --export --models KbS_LV95_V1_4;KbS_Basis_V1_4 --dbfile kbs.gpkg kbs.xtf

passes without errors, but when reimporting the XTF-file the features only appear with point geometries. I suspect that I am missing some export option that states the data consists also of polygon geometries.
For creating the scheme I use:

java -jar ili2gpkg.jar --schemaimport --dbfile schema.gpkg --models KbS_LV95_V1_4;KbS_Basis_V1_4

I am using version 4.3.0 (2019-09-20).


texts: "Interpolation Toolset" is inactive in ArcGIS for Desktop 10.1

texts:
I have ArcGIS for Desktop 10.1 installed on my laptop and "Interpolation toolset" and "Raster Interpolation toolset" is Inactive.

How can I get them work?


texts: query.where statement using date range in JS API fails with IE 11 but works fine with Chrome and Firefox

texts:
My aim to display point data from a feature layer on ArcGIS for JS map interface. I am successful doing so in Chrome and Firefox, however when I include date in a where condition of a query execute it through QueryTask.execute function I am running into issues with IE 11. It works fine when hard coded but not when I read start and end date from a control.

More information:

queryTask.execute throws up the error “Error: Failed to execute query.” on Internet Explorer 11 and EDGE but works fine on Chrome and Firefox browsers. The where condition works fine if I do not give date condition at all in the query for IE11, firefox, chrome, etc.

-----------Case 1:--------------------------------

Case 1: Works fine on IE and other browsers when my where condition is hard coded

query1.where="outlier = 0 and timestamp >= '1/24/2016 5:00:00 AM' and timestamp < '1/25/2016 5:00:00 AM'";
queryTask1.execute(query1, function(results1){.............}, function(error){........});

------------------------Case 2:------------------------

Case 2: is not working on IE. However, works well on Chrome and Friefox browsers.

var WhereCoundition = "outlier = 0 and timestamp >= '1/24/2016 5:00:00 AM' and timestamp < '1/25/2016 5:00:00 AM'";
query1.where = WhereCoundition;
queryTask1.execute(query1, function(results1){.............}, function(error){........});


texts: Merging two rasters, obtain only highest values, overwrite raster in ArcGIS Desktop?

texts:
In order to run a hydrologigal analysis I'd like to to merge two raster datasets (DEM and  buildings) to get as result one raster with only the highest values (elevations) of both. In other words: To overlay the DEM with the buildings and to use the buildings-values whenever possible.

I'm working with ArcGIS 10.4 (ArcMap/ArcScene) for this.

I've got the following input data:

DEM as raster
builidungs as raster (with Null-values outside the buildings)

I was wondering if its possible to use an IF-statement in the Map Algebra Tool, but didn't get the right syntax or maybe there is a better way to do it.
The If-statement could be similar to this:
If ("buildings1" IS Null, "DEM1")


texts: Can ST_Segmentize create equal length segments?

texts:
SELECT ST_AsText(ST_Segmentize(ST_GeomFromText('LINESTRING(0 0, 10 0)'), 3)); gives LINESTRING(0 0,3 0,6 0,9 0,10 0) but the expected result is LINESTRING(0 0,2.5 0,5 0,7.5 0,10 0). How can this be done?

See for example the documentation for qgis:densifygeometriesgivenaninterval:

  "E.g. specifying a distance 3 would cause the segment [0 0] -> [10 0]
  to be converted to [0 0] -> [2.5 0] -> [5 0] -> [7.5 0] -> [10 0],
  since 3 extra nodes are required on the segment and spacing these at
  2.5 increments allows them to be evenly spaced over the segment."


texts: How can I identify least cloudy image in Google Earth engine?

texts:
I need to Load Landsat 8 TOA images, get the least cloudy 2015 image.
Search engine shows several images.

How can I identify least cloudy image?


texts:
I am attempting to create a tool through ModelBuilder that creates address points from parcel polygons. I have added the following tools to the Model. Feature to Point -> Append the point data to our address point feature class. I have been successful with the tool.

I wanted to know if there is a way to move the points after they are appended by let's say 000000.1 of an inch within the tool. I need the points to be displaced by a small amount so that the remaining attribute fields are picked up by ArcObjects (Zip Code & City Boundary). I know the other option is to use the "move command" within the editor toolbar. To have this run in one Python tool would be the ideal result. I have attached the exported Python Code from ModelBuilder, as well as the screenshot of the current model.

# Import arcpy module
import arcpy

# Local variables:
FIREADM_Parcels_MCT = "PARCELS MCT UPLOAD"
Points = "%scratchGDB%\\Points"
FIREADM_AddrPoints = Points
FIREADM_AddrPoints__3_ = "FIREADM.AddrPoints"

# Process: Feature To Point
arcpy.FeatureToPoint_management(FIREADM_Parcels_MCT, Points, "CENTROID")

# Process: Append
arcpy.Append_management("%scratchGDB%\\Points", FIREADM_AddrPoints__3_,
"NO_TEST", "STREET_NUMBER \"STREET_NUMBER\")

 def XYsetVALUE( shape, X_value, Y_value):
   myMoveX = 0.001
   myMoveY = 0.001
   point = shape.getPart(0)
   point.X = X_value + myMoveX
   point.Y = Y_value + myMoveY
   return point


texts: Is there a script to Python for Click to toggle selectable?

texts:
I cannot write a script using ArcPy (Click to toggle selectable).

To make a layer unavailable for sampling!

I just need to create a tool to automate several processes, one of them is the ability to automatically turn off the sample multiple layers on the map! To do this manually! It turns out I have a database with all the layers, I need to when starting a new mxd, clicking on a script or a tool to fully customize my card as I need))

This script does not work, maybe someone else will?

import arcpy
names = "[R_01_001]"
mxd = arcpy.mapping.MapDocument("current")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
layers = arcpy.mapping.ListLayers(mxd)
for layer in layers:
    if layer.name in names:
        unselectable_objs=True
        print unselectable_objs
arcpy.RefreshTOC()
arcpy.RefreshActiveView()

This script does not work

import bpyunselectable_objs = []
for obj in bpy.data.objects:
if obj.hide_select == true: unselectable_objs.append(obj)
obj.hide_select = falseprint(unselectable_objs)


texts: Joining Attributes from Layer Limitation in QGIS 3.2?

texts: I want to join two txt-tables as a layer in QGIS 3.2. In earlier versions (QGIS 2.12) works to join one of the tables more times to the other table. Now in version 3.2, I can join only one time. The table that should be joined second-time is not shown in the drop-down menu to select.


texts: How to create a distance from road raster without using Euclidean Distance tool in ArcMap 10.3.1

texts:
I have created a set of polyline roads as a shapefile, and I would like to create a raster file which shows distance from the nearest road. However, I do not have access to the Euclidean Distance tool so cannot use this. Are there any other suggestions?

I tried buffering the lines then converting this to a raster but this did not work.


texts: Check if a GeoPosition is within a city in OSM

texts:
i am working with GPS Tracks of cars and i would like to find out whether a certain gps point (given in lat/lon) is within a city (and its boundaries/suburbs/etc.) or if it is on the "countryside".
With "countryside" i mean everything that is not defined as "city". Therefore i consider towns, villages, etc. as "countryside".
Basically i want to filter out city-traffic from my data.

Is there a simple way to query this via OpenStreetMaps?

I've already tried to find some tags or attributes on OSM that would indicate whether a gps point is within a city or not. The only relevant attribute i found would be in the "Enclosing Features" tag:
If some point is NOT within a city (and therefore on the countryside) there always seems to be a tag called "City Boundary".
This tag does not exist if the point is within a city.

Example:
Some point (in Berlin) that I would classify as "City":
http://www.openstreetmap.org/query?lat=52.4662&lon=13.4824

Some point (outside of Berlin) that I would classify as "Countryside":
http://www.openstreetmap.org/query?lat=52.3362&lon=13.9618


texts: Taking selected value from one table and populating field in feature class using ModelBuilder?

texts:
I'm still very new to ModelBuilder and I'm trying to work my way through it. I am trying to take a value from a table and basically paste it into a field of a feature class using ModelBuilder. I have a polygon layer that has building footprints and another polygon layer that is squares for flooding depth. These two layers can not be joined.

I currently have the model set up to select 1 building at a time and then use that building to select by location via intersect on the flood layer. This then selects anywhere between 1 and 6 flood squares. I then have to take the average of the flood depth from those squares and add it to an existing field in the building layer. This is the part I get stuck on. I can run a summary statistics to get the average of the selected flood square which creates another table.

How would I get that value into a field in the building footprint layer?

If the calculate field could accept multiple inputs this would be simple but I am having a hard time taking the data from one table and filling a field in a different shapefile.


texts: It seems there is a bug in St_Clip function of PostGIS

texts:
The code I am using for clip feature is :

raster ST_Clip(raster rast, geometry geom, double precision[] nodataval=NULL, boolean crop=true);

I tried clip a shape depends on the polygon from the Raster Data. But the PostGIS would clip the raster depend on the Bounding Box of the polygon, and put the Non-data value as 0. Even I tried to modify that with code:

SELECT ST_Clip(rast,geometry, ARRAY[-9999],true) FROM dtm, polygons WHERE polygons.polygonid = 1;

After I dump the data as a ASC-II file from the PostGIS, using gdal_translate with another option -a_nodata -9999, I find out it only changes the Header of the ASC-II file not the exact bounding box value.

Is there anyway to change the value from 0 to -9999, or someway to clip the raster depending on the polygon rather than the Bounding Box of the polygon.

THanks.


texts: Manual (or examples) of r.ros, r.spreadpath and r.spread in GRASS

texts:
you can report me references or applications for use of the following modules (GRASS):
- r.ros,
- r.spreadpath,
- r.spread.


texts: Using the Proj4js library to convert from Google Maps to projected values

texts:
I'm using the proj4s javascript library (based on advice from this question I asked previously) to convert lat/long values on my webpage. This page receives the lat/long values from google maps and converts them to the British National Grid (ORD SURV GB), which wikipedia tells me has the following code: EPSG:27700. I want to calculate the maximum and minimum values of both the latitude and longitude, so I can query a database.

My code is as follows:

var source = new Proj4js.Proj('WGS84');
var dest = new Proj4js.Proj('EPSG:27700');

bounds = map.getBounds()  // Get bounds of map after user pans around
zoom = map.getZoom()      // Get zoom level

xMax = bounds.getNorthEast().lat()  //get max and min values for lat/long
yMax = bounds.getNorthEast().lng()
xMin = bounds.getSouthWest().lat()
yMin = bounds.getSouthWest().lng()

alert('xMin: ' + xMin + '\nxMax: ' + xMax)
alert('yMin: ' + yMin + '\nyMax: ' + yMax)

// transforming point coordinates
var latXNew = new Proj4js.Point( xMin, xMax );
var latYNew = new Proj4js.Point( yMin, yMax );

Proj4js.transform(source, dest, latXNew);
Proj4js.transform(source, dest, latYNew);

However, I'm not getting the correct results from my webpage. I've loaded in a kml file that has a rectangle with boundaries that I know, and the values don't match up. They appear to be off by at least a factor of approximately 8.

NE Lat: 52.25635981528854
NE Lng: 1.3534606328125243

Expected Lat: 500000.0
Expected Lng: 150000.0

Returned Lat: -5527596.989640327 (incorrect)
Returned Lng: 622727.1071581601 (incorrect)

What am I doing incorrectly?


texts: Printing Labels of Unique Value Classification Using ArcPy?

texts:
I have a layer where I have symbolized the features on unique values. I have also grouped several features together. EX) Layer: Land Use; Value Field: Type of Land Use; Unique Values of different park types have been grouped together and there 'Label' changed to "Parks"; All other unique values left the same.

Using ArcPy 10.2, how can I print the "Label" (not the unique value itself) of the Unique Value categorization? I am looking for something where I can use arcpy.mappings?


texts:
In Iso4app QGIS plugin, in the isochrone parameter window, how speed type and speed limit work?
do they work based on the maximum allowed speed of the roads? or if we determine a speed level more than maximum allowed speed of the road in the speed limit field it would create isochrones based on the entered speed?


texts: Creating polygon out of list of points QgsPointXY

texts:
I am attempting to create a simple polygon on QGIS 3.14 using PyQGIS by adapting the code from How to draw polygons from the python console?.
I am running into an error with the setGeometry line where the function fromPolylineXY appears to need a QgsPointXY object but the pts object is a list.
Code so far is listed below.
tmp = QgsVectorLayer('Polygon?crs=epsg:29194', '200905_Bdy',"org") #Layer for polygon
prv = tmp.dataProvider()                               #Data object for polygon layer
ply01 = QgsFeature()                                   #Object for polygon
# Object with polygon verticies
pts = [QgsPointXY(396100,8969000),QgsPointXY(396100,8973900),QgsPointXY(397900,8973900),QgsPointXY(397900,8969000)]
ply01.setGeometry(QgsGeometry.fromPolylineXY([pts]))


texts: Geoserver WFS GetFeature query that matches entire word in string

texts:
I am trying to write a WFS getfeature query for some GeoServer JSON data that returns features that match a part of a string only if the entire word matches. The word can be in the start, middle, or end of the string. I started by using cql_filter, but am returning records that don't match a whole word.

Say I have two features with the following values for AddLocID:

id, AddLocID
1,  216
2,  216 72216
3,  72216

I want to return features 1 & 2 and not feature 3

Using cql_filter=AddLocID=216 returns only feature 1 and does not match feature 2

https://bloomington.in.gov/geoserver/publicgis/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=publicgis:HistoricSitesandStructuresSurvey2018&cql_filter=AddLocID=216&outputFormat=application%2Fjson

Using cql_filter=AddLocID Like 216 returns all 3 features but I only want 1 and 2

https://bloomington.in.gov/geoserver/publicgis/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=publicgis:HistoricSitesandStructuresSurvey2018&cql_filter=AddLocIDlike%27%216%%27&outputFormat=application%2Fjson

How can a write my query so that records 1 and 2 are return that match 216 as a whole word in the string value?


texts: How to use Advanced Symbology in ArcGIS for Android?

texts:
I try to use the advanced symbology in order to draw nato symbols in JAVA
I need little help about it,

What the parameters _WKID mean like in this command : lmessage.setProperty("_WKID", "3857");
What the parameter UniqueDesignation mean?
THE MOST IMPORTENT QUESTION: How can send parameters to symboles llike the symbole "Target" (TACGRP.FSUPP.PNT.TGT.PTGT) (2.X.4.1.1.1) with sic G*FPPTS---**X
parameters like "Target Alt" , "Description" and "Name"?


texts: Script or model to create random points having multiple minimum distances within a polygon

texts:
Given:  A single part polygon, total number of points to place inside that polygon, and table specifying the size of each subset (number of subsets n >= 1) of that total number of points and each subset's minimum spacing.

Find:  A (semi-)random distribution of points inside the polygon such that no point is closer than its minimum spacing to any other point in the polygon.

For example, you are given a polygon and 100 points.  Of the set of 100 points there are three disjoint subsets--A, B and C.  Subset A contains 50 points having a minimum spacing of 7 m, subset B contains 25 points having a minimum spacing of 4 m, and subset C contains 25 points having a minimum spacing of 2 m.  Distribute the points such that no point in subset A is within 7 m of any of the remaining 99 points, no point in subset B is within 4 m, no point in subset C is within 2 m.

One possible method to solve this would be to start with subset A, distribute the points randomly using Create Random Points, create a 7 m buffer around the points, erase those buffers from the original polygon and use the erased polygon to randomly distribute the points in subset B. . .repeat through all subsets.

I'm looking for a way to implement this (or a superior) solution using scripting or model builder.  Ideally, the user would input the number of points in each subset, the minimum spacing for each subset, and the polygon feature class into which the points would be placed.


texts:
I want to create a multimodal network dataset for the city of Montreal. I have roads and transit lines features for Quebec province (where Montreal is). Because I need the features limited to Montreal only, I clipped them using Clip tool in geoprocessing menu with Montreal census tracts layer as the clipping feature. Now, the problem is that Montreal has water bodies which are not a part of census tract layer. The areas of water bodies are therefore blank spaces between census tracts. The clipping procedure cuts the roads and transit lines where they pass over the water bodies. This, obviously, can't be used for creating an accurate network datset. Which tool can I use to clip the road and transit lines feature classes to city of montreal census tract layer only, without cutting them?
Alternatively, is it possible that I create network datset for the whole province and then clip it?


texts: Finding nearest 10 observations that occurred 2 years before original observation using ArcGIS Pro ModelBuilder?

texts:
I have a point sample of house sales (~42k points over 9 year period) with a date field(integer), and for each point I want to find the 10 nearest houses within a year of sale and also which were sold before the house in question.

The plan is to use the iterate feature in ModelBuilder to select each point, and also the sample of near features, run the generate near table for the 10 closest houses, append to a table and label using the value output field from the iterator, and loop again.

How do I generate the near features sample for each individual house?

The table is ordered by objectid=1 is the most recent sale and the highest objectid is the oldest sale.


texts: Configuring manual python library installations

texts:
Windows 7.  My python path reads:

C:\Python27\ArcGISx6410.4;C:\Python27\ArcGIS10.4;C:\Anaconda;C:\Anaconda\Scripts

The end of my user path reads exactly the same thing.

Yet, when I try to install libraries - for example, the requests library:

python setup.py install

I get:

C:\Python27\ArcGISx6410.2 not found

I don't understand. I uninstalled ArcGIS 10.2, have removed all current user registry keys (I don't have full admin privileges), and my path variables don't reference the 10.2 installation at all.

When I bring up the python console in ArcMap and then type "import requests" it says "No module named requests."

What I need is to configure my command line installs to by default go into C:\Python27\ArcGISx6410.4 for accessibility within Spyder,

and use the --prefix install option to send them to C:\Python27\ArcGIS10.4 to access the same libraries from within ArcMap, which confusingly only uses 32 bit in foreground.


texts: Performing interpolation with QGIS in Python application?

texts:
I have to make an interpolation map for river water quality, but it must be made in the Python programming language.

Can I use QGIS in a Python application?

Not Python in QGIS but QGIS in Python.


texts: Opening embedded canvas in new window with all layers

texts:
I am creating a standalone custom PyQGIS application and have some questions regarding canvas.
I have read the cookbook:
https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/canvas.html
They mentioned that you can open the canvas with
canvas = QgsMapCanvas()
canvas.show()

However, I have already embedded the canvas into my GUI. I wish to have the embedded canvas with all my layers in it pop up as a separate window when I click a button. I have tried this:
self.canvas = QgsMapCanvas(self.centralwidget)
self.canvas.show()

But it doesn't work.
QtDesigner test Code:
from PyQt5 import QtCore, QtGui, QtWidgets

class Ui_MainWindow(object):
    def setupUi(self, MainWindow):
        MainWindow.setObjectName("MainWindow")
        MainWindow.resize(800, 600)
        self.centralwidget = QtWidgets.QWidget(MainWindow)
        self.centralwidget.setObjectName("centralwidget")
        self.horizontalLayout = QtWidgets.QHBoxLayout(self.centralwidget)
        self.horizontalLayout.setObjectName("horizontalLayout")
        self.pushButton = QtWidgets.QPushButton(self.centralwidget)
        self.pushButton.setObjectName("pushButton")
        self.horizontalLayout.addWidget(self.pushButton)
        self.widget = QgsMapCanvas(self.centralwidget)
        self.widget.setObjectName("widget")
        self.horizontalLayout.addWidget(self.widget)
        MainWindow.setCentralWidget(self.centralwidget)
        self.menubar = QtWidgets.QMenuBar(MainWindow)
        self.menubar.setGeometry(QtCore.QRect(0, 0, 800, 26))
        self.menubar.setObjectName("menubar")
        MainWindow.setMenuBar(self.menubar)
        self.statusbar = QtWidgets.QStatusBar(MainWindow)
        self.statusbar.setObjectName("statusbar")
        MainWindow.setStatusBar(self.statusbar)

        self.retranslateUi(MainWindow)
        QtCore.QMetaObject.connectSlotsByName(MainWindow)

    def retranslateUi(self, MainWindow):
        _translate = QtCore.QCoreApplication.translate
        MainWindow.setWindowTitle(_translate("MainWindow", "MainWindow"))
        self.pushButton.setText(_translate("MainWindow", "Fullscreen"))
from qgis.gui import QgsMapCanvas

if __name__ == "__main__":
    import sys
    app = QtWidgets.QApplication(sys.argv)
    MainWindow = QtWidgets.QMainWindow()
    ui = Ui_MainWindow()
    ui.setupUi(MainWindow)
    MainWindow.show()
    sys.exit(app.exec_())

Main code:
class Main(QtWidgets.QMainWindow, Ui_MainWindow):
    def __init__(self, parent=None):
        super(Main, self).__init__(parent)
        self.setupUi(self)

        self.pushbutton.clicked.connect(self.makeFullscreen)

    def makeFullscreen(self):
        self.widget.show()

The output should be a canvas in a new window after clicking the button.


texts: CORS error when try to display layer from remote geoserver

texts:
I am following this tutorial in order to add vector tiles using geoserver and OL3.
I am using some layer from a remote geoserver which I have access:

var layer = 'myworkspace:arm_trs_roads_osm';
var projection_epsg_no = '4326';
  var vectorLayer = new ol.layer.VectorTile({
  style:style_simple,
  source: new ol.source.VectorTile({
    tilePixelRatio: 1, // oversampling when > 1
    tileGrid: ol.tilegrid.createXYZ({maxZoom: 19}),
    format: new ol.format.MVT(),
    url: 'http://ogcserver.gis.wfp.org/geoserver/gwc/service/tms/1.0.0/' + layer +
        '@EPSG%3A'+projection_epsg_no+'@pbf/{z}/{x}/{-y}.pbf',
    crossOrigin: null
  })
})

When I try to display the layer on the map I get the well known CORS error:

Failed to load http://ogcserver.gis.wfp.org/geoserver/gwc/service/tms/1.0.0/my_workspace:arm_trs_roads_osm@EPSG%3A4326@pbf/3/0/4.pbf: No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'http://vagrant-test.local.com' is therefore not allowed access. The response had HTTP status code 404.

I am not sure how to solve this issue. I found this question which says to enable CORS in the tomcat level. But it also says that this feature is available after tomcat 7.0.41 (I have 7.0.39; not possible to upgrade currently).

I have also enabled JSONP in my geoserver config files. But still I get the same error. Which other alternative do I have?


texts: How to use SAGA_CMD for batch processing?

texts:
I am trying to use SAGA to triangulate a set of .las files I have. I have tiled them using LAStools, since the single file was too big for my computer to handle. I wonder, can I use the saga_cmd grid_gridding "Triangulation" for multiple files, or do I have to process each individually? And I want tif files based on the Z-value, where do I specify the Z value? Is it the -FIELD Attribute?

I have found very little information about the cli scripting of SAGA, are there any book available on the subject?


texts:
I have a geojson with height informations in it. I tried adding the geojson to the scene using

var viewer = new Cesium.Viewer('cesiumContainer', {
sceneMode : Cesium.SceneMode.SCENE3D,
timeline : false,
animation : false
});

var dataSource = Cesium.GeoJsonDataSource.load('../../SampleData/Buildings_3D_Elevation.geojson');
viewer.dataSources.add(dataSource);
viewer.zoomTo(dataSource);

I can view the data as below. What should be done to use the height information and view as 3D?


texts: Adding a Layer Control Legend based on geoJson properties

texts:
I have a Leaflet map app in which I am having difficulty getting the L.Control to function properly. The app pulls points based on data from a geoJson file. The geoJson data has a property "Type" (feature.properties.Type). The types are restaurants, bars, music, cafe, and culture. I want a Layer control legend that will toggle the points on and off the map depending on whether the use has them clicked in a checkbox. For example, if the user has bars clicked, but not culture, all the points with Type == 'bar' will show on the map but not those with 'culture'. This is how the geoJson is set up; mind that I already have variable to control what the icons look like that I am not including:

var points = L.geoJson(dtr_points, {
onEachFeature: function (feature, layer) {
    var link_url = "" + feature.properties.Name + ""
    console.log(link_url);
    layer.bindPopup(link_url);

     if(feature.properties.Type == "beer") {
        layer.setIcon(beer);
     };
     if(feature.properties.Type == "bar") {
        layer.setIcon(bar);
     };
     if(feature.properties.Type == "restaurant") {
        layer.setIcon(restaurant);
     };
     if(feature.properties.Type == "music") {
        layer.setIcon(music);
     };
     if(feature.properties.Type == "culture") {
        layer.setIcon(culture);
     };
     if(feature.properties.Type == "cafe") {
        layer.setIcon(cafe);
     };
     if(feature.properties.Type == "occult shop") {
        layer.setIcon(occult);
     };
  }
}).addTo(map);

L.control.zoom({
  position:'bottomright'
}).addTo(map);

L.control.layers(points).addTo(map);

Notice that the last line is where I am trying to add the Layer Control. It does not show up on the map currently. Can you suggest how to fix this? I think I need to add a layer group, but the Leaflet documentation doesn't explain how to do this for geoJson data.


texts: Calculating distance between multiple lines of points in QGIS?

texts:
I'm looking to calculate a distance from a sensor reading to that of the nearest sensor reading on a different sensor line, all held within the same points file. Ideally using QGIS.

The image below is an attempt to picture the problem: the red dots are our points. Black lines are to show the points are associated with that line. The grey distances are what we want to associate with the points.

I have looked at nearest neighbour and hub distances but struggling to work out a workflow.


texts:
I installed ArcGIS for Desktop 10.4.1, ArcGIS for Server 10.4.1 and Web Adaptor, now I need to install Portal for ArcGIS, I looked for the trial version in the ESRI's website but in vain.

Is there any link I can access to get the Portal?


texts: Python: Break linestring based on condition

texts:
I have a geopandas dataframe of a bunch of linestrings that have some data associated with each vertex/point:

Point_x = (Lat, Lon, Time, ID, Data1, Data2, Data3)

The points are converted to linestrings based on ID and ordered by Time.

I want to break the linestrings where at the point where some condition is met. Right now that's when the distance between Points is greater than some value. In the future it could be where a function of the Data fields is some value. For instance, split a linestring when Speed crosses 5 kph.

The current problem is that some of the tracks are formed from points that have duplicate ID's so the linestring jumps back and forth over huge distances and I want a threshold to break these lines.

Any ideas on the correct way to structure this or libraries/methods that might be of use?

The dataframe has over 150k tracks with many points per track in it so efficiency would be nice.

Here's an example of the tracks DF:

ID         geometry
204235000  LINESTRING (37.62001 -28.99535, 37.62015 -28.9...
205400000  LINESTRING (3.807816666666666 -18.083181666666...
207138000  LINESTRING (22.73206 -34.97915833333333, 22.73...
209016000  LINESTRING (8.447673333333331 -23.522783333333...

Here's an sample from the points DF. There are 18 columns including Datetime, Point(Lon, Lat), Speed, Size etc etc:

Index           Heading   Latitude  Longitude       ID
20              92.8 -35.946802  13.089695  210725000
21              93.5 -35.946912  13.091808  210725000
22              95.4 -35.965520  13.497698  210725000
23              94.7 -35.965803  13.501898  210725000
24              94.9 -35.965987  13.504573  210725000

EDIT: Tried to be a little clearer.


texts: QGIS: How to create attribute following polygons' by object-to-object relationships?

texts:
I'm trying to map the sub-basins within North America's Great Basin.  I'm using the National Hydrography Dataset, 12-digit watershed boundary data.  Each polygon has an ID, also either the ID of the downstream polygon, or a code denoting it a terminal basin.
In QGIS, how do I create an attribute populated by the ID of the eventual terminal polygon of each polygon?
For example, consider the following attribute table.  The ID and TO_ID attributes exist; I want to generate the TERMINAL_ID attribute.  The value T designates a terminal polygon.

ID (exists)
TO_ID (exists)
TERMINAL_ID (generate)

1
2
5

2
3
5

3
5
5

4
5
5

5
T
5

6
8
8

7
8
8

8
T
8

9
T
9


texts: Use icons for polygons based on Boolean values from attribute table as secondary symbology

texts:
I have a multypolygon-shapefile of government districts in Germany with two columns, "population density" and "yesno":
The first has integer the second one Boolean values.
I have already colored the map using the population density values (Symbology - graduated).
Now I'd like to use the second column (yesno) for symbology as well.
I want to put one of two symbols (one for true and one for false) on every government district.
I know how to use icons for points, but not for polygons.
And I know how to color based on one column, but not two for the same layer.

How can I use an icon for a polygon without losing the earlier classification?


texts: Iterate through raster in R and apply conditional statements

texts:
I have a raster data set in R with numeric values ranging from 0 to 61.5, including some random NAs throughout. I would provide a reproducible data set, but I don't know how to create a raster with random NAs and I cannot share my actual data online. I simply need to iterate through each raster value and apply the following logic:

> for (i in 1:length(raster){
>  if (value(i) >= 5.1) {value(i) = 1}
>  else if (value(i) < 5.1) {value(i) = value(i)/5.1}
>  else (is.na(value(i)) == TRUE) {value(i) = 0.000001}
> }

I often do loop on data frames in R but I'm new to the raster format and can't quite get this up and running. I'm not sure how to access each raster value in a for loop and apply my conditional statements. If something isn't clear let me know in the comments and I will fix my question as needed.


texts: alternative to geograpy package for python (Pulling locations from text)

texts:
I am currently exploring using the geograpy package for python to pull locations from text. It uses nltk, among other things to pull locations from URL's and text.

The only problem is it doesn't support two word places, like "South Carolina".

Is there a good python package/gazetteer to try to pull these? I am trying to feed it into a twitter feed mapper.


texts: QGIS to have infinite map canvas in East-West direction

texts:
can I setup QGIS in a way to scroll infinitely in east-west direction?
To achieve this, QGIS would need to repeat the underlying chart (XYZ Tiles or WMS) once the boundary is reached.


texts: Ordered weighted averaging error: can only concatenate list (not "str") to list

texts:
I am trying to perform an ordered weighted average on two rasters in QGIS 3.0.3 (as a trial run for doing it on more than two rasters). But I keep receiving the following error code:

Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python/plugins\processing\algs\saga\SagaAlgorithm.py", line 275, in processAlgorithm
s = values[i] + '\t' + values[i + 1] + '\t' + values[i + 2] + '\n'
TypeError: can only concatenate list (not "str") to list

Execution failed after 0.05 seconds

I googled the error message but my programming knowledge is limited.
I first thought that it could be because the rasters were of different types so I converted both to UInt16 (Raster->Conversion->Translate) but I am still getting the same error. I am not sure what to try now to fix the problem.

Here is some background information on the rasters and the project I am working on:

I am trying to produce an overall risk map considering the relative contribution of various factors (e.g. rivers, roads).
Each raster contains buffer zones around one factor (e.g. roads).
Each buffer zone is assigned an integer value representing the level of risk (0= no data, 1= no risk, ..., 7=high risk).
For this trial run I was attempting to weight two layers at 60, 40.
One layer had integers from 0-5, the other had integers from 1-7.
The rasters are single band, and tif format.


texts: Measuring distance between related points using ArcGIS Desktop?

texts:
How can I generate a table with measured distances between two related point layers?

Every point from A layer has at least 1 assigned point in B layer and the attribute that defines relation between them is saved in separate table (A layer->table-> B layer). I have made relations in arcGIS to define pairs of points. Now I need to know what is the distance between them.


texts: QGIS - is it possible to select the columns shown with "Identify features"?

texts:
In QGIS, I have created a join to an attribute table with some 60 entries. I am only interested in seeing 3 of the columns, specifically address data.
I can see no way of excluding the unwanted columns from the “identify features” tool.
Can anyone help?


texts: Are there working build instructions or prebuilt binaries of gdal stable >= 2.3.0 with pdfium available for linux x64?

texts:
I have tried to manually download gdal latest (2.3.0 at time of this post) so that I can compile it on linux x64 (ubuntu 16.04) with the --with-pdfium flag enabled to test out GeoPDF to MapBox Vector Tile conversion. After much tinkering I seem to have successfully compiled pdfium with https://github.com/klokantech/pdfium/blob/build/build-lin.sh, but when I go to compile gdal it errors out with undefined reference to CFDF_Document::CreateNewDoc().

Has anyone successfully built gdal 2.3.0 or above with the --with-pdfium on linux? Can you link to working instructions? Are there precompiled binaries available?

Edit: I can't get this working so I'm just going to use podofo which is installable from apt-get


texts: How to color polygons based on attribute in QGIS 3.8?

texts:
I want to color polygons (fill color) based on an integer attribute (named "fold", ranging 1-5). It was so simple in older versions of QGIS (2.x), but now I installed version 3.8 and I am clueless (see the screenshot). I went to Layer Properties as usual, looking for the Style tab, which is no longer there, so I ended up in similar Symbology tab, but I don't see a way how to do it here.

PS: giving up, had to go back to the simplicity of good old QGIS 2.8


texts: Estimating cline width at multiple points along a trait gradient

texts:
In the map below I have 66 (green) points along the midpoint (0.5) for an allele frequency cline across the region (see map). I created the points by placing a point every 25 km along the 0.5 allele frequency contour. The thick red lines represent the 0.22 and 0.78 allele frequency contours. I would like to estimate the width of the cline at each of the 66 points. The "estimated" width should be based on the shortest distance between the two contours when the line runs through the point. This would allow me to compare the cline width at 25 km increments.

The result would look something like the map below. Note that I hand-drew only of 3 of the 66 desired lines and my lines are probably not the shortest distance between the two contours for each of the 3 points. Each of the line must be a straight and unbroken (i.e. no bend at the point). I would also get the total distance of the line, plus the distance for each segment (0.22 to 0.5 and 0.5 to 0.78).

This is a map created using the "Near" function with my points as the "input feature" and contours as the "near feature". I created the map by running the near function for each contour separately. As you can see the map is not quite right. Any other suggestions?


texts: Why does an 32-bit FGDB raster have a smaller file size than 8-bit?

texts:
Using Arc 10.0, I have a single band, 32-bit signed raster as the only item in a file geodatabase.  Using ArcCatalog, I copy that raster to a new geodatabase with an 8-bit unsigned pixel type.  As expected, Raster Properties shows that the 8-bit raster is smaller (438 MB) than the 32-bit raster (1.71 GB).  However, when I measure the size of the two gdb folders using Windows File Explorer properties, the 32-bit folder is smaller(!) than the 8-bit folder (64 vs 83 MB, respectively).  Both gdb folders contain 66 files.

Both rasters are 32,225 X 13,417, single band, 30 meter cell size, FGDBR format, continuous Source Type, no colormap or pyramids, with LZ77 compression.  The only difference is the pixel depth and type: 32-bit signed integer vs 8-bit unsigned integer.

Why is the 32-bit gdb folder smaller than the 8-bit folder, as measured by Windows File Explorer?
Why does Arc raster properties report such wildly inflated size (438 and 1,710 MB) compared to Window File Explorer (83 MB and 64 MB)?

In ArcMap, the 32-bit image takes about 9 seconds to display, while the 8-bit image loads almost instantly.


texts:
I am following the book "Geoprocessing with Python" and I have encountered an issue with my code:
from osgeo import ogr

shp_ds = ogr.Open('path/to/data')

volcano_lyr = shp_ds.GetLayer('us_volcanos_albers')

cities_lyr = shp_ds.GetLayer('cities_albers')

#How far is a particular city from a volcano?
volcano_lyr.SetAttributeFilter("NAME = 'Rainier'")
feat = volcano_lyr.GetNextFeature()
rainier = feat.geometry().Clone()

cities_lyr.SetAttributeFilter("NAME = 'Seattle'")
feat = cities_lyr.GetNextFeature()
seattle = feat.geometry().Clone()

#Distance function:
meters = round(rainier.Distance(seattle))
miles = meters / 1600
print('{} meters ({} miles)'.format(meters, miles))

I get this error:
  File "[redacted]", line 136, in <module>
    seattle = feat.geometry().Clone()

AttributeError: 'NoneType' object has no attribute 'geometry'

It opens perfectly fine in QGIS and the columns are named the same. Additionally, I managed to analyse how many cities there are in the city layer (redacted here), so this file should definitely have geometry.
Would anyone be able to assist?


texts: Google Earth Engine export CRS option definition

texts:
I have been using the Global Surface Water dataset for some analysis, but export the data to my local machine before using it. The dataset is unprojected on a WGS84 datum (EPSG 4326), but I want to analyze it in the native Landsat CRS, i.e. the respective UTM zone of each study site. When I use export, I'm simply entering the target CRS in the "Export" call. This is causing a couple of issues. When I compare the exported scene with the scene available in GSW, it looks like averaging between landsat 7-derived no data bands and existing water bodies (like channels) leads to spurious land pixels. In addition I notice significant warping of objects like lakes. An example comparison is below. The coordinates of this location are 69.601983 deg N, 161.742794 deg E.

This leads to the two questions I have:

1) Is specifying the CRS a reprojection from EPSG 4326 to EPSG X or redefining the CRS of the exported map (incorrectly)?

2) If this is a reprojection what is the reprojection method? The GSW data is categorical but takes on integers 0 to 2 (0 = No data, 1 = Land, 2 = Water), so a nearest neighbor resampling won't cause the above issue, but averaging or interpolation will.

The script I use to export over a given region is shown here. I should note that this region is on the edge of UTM 57 N and UTM 58 N so part of the warping may be because of that (but hopefully not). I have compared using 58 N in the export and get similar issues. I know that one solution to this may be to mark the 0 values with NA (or whatever JS equivalent) before exporting, but the significant warping is still bothersome.

UPDATE

I found that the reprojection in GEE appears significantly off from what I thought I was downloading. The script below was used to download data in EPSG:32657 (UTM57N). Alternatively, I downloaded data specifying no CRS and a scale of 30. This exports an image in EPSG:4326 with res = .000269 deg. I tested reprojecting to 15 m and 30 m resolution UTM57N scenes and found much better correspondence (visually). Four scenes are shown below which highlight the issue with the GEE export. I understand the issue is when I pass "30" as an argument to scale in export and the crs as the UTM zone, so to tie it back to Question 1: this is probably a reprojection but what is the new scale argument doing?

/// Visual parameter constants
// Black signifies no data 0, grey is land 1, white is water
var VIS_WATER_MASK = {
  min: 0.0,
  max: 2.0,
  palette: ['black', 'grey', 'white']
};

///RECURRENCE
//var date = "recurrence_07";

//var grab = ee.Image('JRC/GSW1_1/MonthlyRecurrence/monthly_recurrence_06');
print(grab)
/// Selection of individual months
//////////////////////////////////////////////// May
 var date = "2000_07";
 print(date);
 var filename = 'JRC/GSW1_1/MonthlyHistory/' + date;
 print(filename);

 var grab = ee.Image(filename);
 print(grab)
  var water1 = grab.select("water");
 Map.addLayer({
   eeObject: water1,
   visParams: VIS_WATER_MASK,
   name: date
   });

Export.image.toDrive({
  image: grab,
  description: date,
  scale: 30,
  region: Kolyma,
  crs: 'EPSG:32657',
  folder: 'Kolyma_GSW_32657',
  maxPixels: 9e8
});

This is a GSW mask in EPSG 4326 (i.e. unprojected) with .000269 deg resolution, exported with scale = 30. This is about 15-m a pixel here (at 69 N).

This is if you use the simple script I put above to export, setting CRS = 32657 (UTM 57N) and scale = 30.

This is when doing projection from 4326 to 32657 at 15 m spatial resolution (automatically set without specifying a resolution) using QGIS instead of GEE:

This is when doing projection from 4326 to 32657 at 30 m spatial resolution (prespecified target resolution) using QGIS instead of GEE:


texts: Display both data in different fields in QGIS attribute table

texts:
In my attribute table I have lat_1, long_1 and lat_2, long_2 in the same row. How can I display both at the same time. Is their any solution for that in QGIS.
lat_1, long_1 are main location
lat_2, long_2 are entry point location
That means in attribute table field names should be unique right. So it will display only x and y. But in my attribute table their are entry points (lat, long) they are also coordinates. So I want the both x & y (main and entry) displayed.


texts: Python script in ArcGIS to add field to a feature class in a workspace

texts:
I'm facing a difficulty to identify a feature by its name in a workspace and addfield to it.

import arcpy
arcpy.env.workspace ='C:\DATA_GEO\Bases de données\Vague 2\BOURRHA'
fcList = arcpy.ListFeatureClasses()
for feature in fcList:
    if feature.rfind("centroid")!=-1: #look for the feature having "centroid" in its name
        arcpy.AddField_management("fc", "ref_ID", "LONG", 9, "", "", "refcode", "NULLABLE", "REQUIRED")

But this is the error Python returns

  Runtime error
  Traceback (most recent call last):
    File "", line 3, in
  TypeError: 'NoneType' object is not iterable

How can I solve this?


texts: Select layer by location and copy features from multiple shapefiles to multiple shapefiles

texts:
I want to apply the reverse of the intersection tool between two shapes several times for different shapefiles. I have 95 point features in folder A, each representing a species location. In folder B, I have 95 polygon buffer areas, each representing potential locations for the 95 species. For instance, I want to pair species 1 in folder A, with species 1 in folder B, species 2 in A with species 2 in B and so on. From there I want to get the points that do not overlap between every pair of features.

In ArcPro I use select layer by location and set intersection as the spatial relationship, then check "invert spatial relationship", finally I generate a copy of the features that do not overlap.

For every species, I tried something like this but still need to define a method of pairing for every species in the two folders:

In selecting features I added all 95 species in a single feature shapefile

In R the code for a sinlge pair of shapefiles is as follows:

set_here()
library(here)
library(sf)
library(dismo)
library(spex)

presence.data <- st_read(here("data/presence_data.shp")) %>%
  st_transform("+proj=utm +zone=17 +datum=WGS84")

hs_spp <- st_read(here("data/hs_spp.shp")) %>%
  st_transform("+proj=utm +zone=17 +datum=WGS84")

presence.buffer <- st_buffer(presence.data, dist = 5000, endCapStyle = 'ROUND')

absence.data <- st_difference(hs_spp, st_union(presence.buffer))

st_write(absence.data, here("/scratch/abs.shp"))

How can I generate an output future shapefile that represents the points not overlapping with a buffer for every 95 species?


texts: Change language of QGIS Browser to English

texts:
I want to change the language of QGIS Browser graphical user interface (GUI) to English. I know how to change the language of QGIS GUI itself from (Settings->Options->Locale), but that does not change the language of QGIS Browser. I also tried to change the language of QGIS Browser from OSGeo4W Shell with the following commands:

qgis-browser --lang=en
qgis-browser --lang en
qgis-browser --lang en_US
qgis-browser.bat --lang en_US
qgis-browser.bat --lang=en_US

But unfortunately none of the above commands succeeded to change the language of QGIS Browser.

Do you know how can I do it?


texts: I'm using QGIS 3.14 on Windows 10.   When I try to use the 'identify features' tool, it clearly highlights the item I'm looking to identify, but no results panel appears.  There is nothing docked and it doesn't seem to matter what item or layer I try to identify, nothing shows up.  I've closed all other windows and tried restarting QGIS, but nothing appears.


texts: Creating heatmap from field/column in attribute table of shapefile in QGIS

texts:
I'm almost completely new to QGIS, so help would be greatly appreciated :).

I have a layer from a shapefile loaded up into QGIS with a few thousand zones/features. There are also multiple fields/columns in the attribute table. I need to create a heatmap (or something like this)

from one of those fields, let's call it X. Field X has a value for each zone/feature.

I already looked around a bit, but if I go to the layer properties --> Symbology, I don't see heatmap in the top dropdown menu.

Can anyone help me with how to create a heatmap?


texts: Constructing a raster with Moran's I close to -1

texts:
I wrote a script that produces three raster objects and computes their local and global values for Moran's I. The rasters exhibit regular, clustered and randomly distributed data. I would have expected the values for Moran's I to approximate -1, 1, and 0 for the three scenarios. However this is not the case for the regular data. Moran's I comes out at around 0 instead of -1, both locally and globally and I don't understand why.
To form this into a question: What do I have to change about the regular grid so it results in a Moran's I approaching -1?
Here's the script that I used to compute and plot the rasters and statistics (cell neighbourhood definition is queen's case):
library (raster)
library (rasterVis)

n <- 11

r_regular <- raster (nrows=n, ncols=n, xmn=0, xmx=1, ymn = 0, ymx = 1)
values (r_regular) <- 0
values (r_regular) [seq (1, n^2, by = 2)] <- 1

m_regular <- MoranLocal (r_regular)
mg_reg <- Moran (r_regular)

r_clust <- r_regular
values (r_clust) <- 0
values (r_clust) [seq (1, (n^2) / 2)] <- 1

m_clust <- MoranLocal (r_clust)
mg_clust <- Moran (r_clust)

r_rand <- r_regular
values (r_rand) <- 0
values (r_rand) [sample.int (n^2, n^2 / 2)] <- 1
m_rand <- MoranLocal (r_rand)
mg_rand <- Moran (r_rand)

ms <- stack (r_regular, r_clust, r_rand, m_regular, m_clust, m_rand)
nms <- c ("Regular data", "Clustered data", "Random data",
                 paste0 ("Regular - Global I: ", mg_reg),
                 paste0 ("Clustered - Global I: ", mg_clust),
                 paste0 ("Random: Global I: ", mg_rand)
)
levelplot (ms, names = nms)


texts: Setting python code to trigger on any update on feature class?

texts: For example, when I have a geometry change on a point, I would like to update the X and Y fields to the new current location, I understand how to make the code, but not where do I set or put it to be triggered/called


texts: Difference between Execute Task and Submit Job?

texts:
When I compare my GPServer with sample service from ESRI. I notice difference.
In esri

Supported Operations:   Execute Task

my GPServer:

Supported Operations:    Submit Job

What is the difference between Task and Job??


texts: WFS layer from GeoServer with OL4 - String index out of range error

texts:
I need to load a WFS layer from GeoServer 2.8 to a web app using OpenLayers 4.6.5 and JQuery 2.2.3 so that I can use Select to highlight a feature on click.
I can preview the line layer as a WMS with no problems in GeoServer, but when I try loading in the app, in the network response gives this error:

  <ows:Exception exceptionCode="NoApplicableCode">
    <ows:ExceptionText>java.lang.RuntimeException:
    java.lang.StringIndexOutOfBoundsException: String index out of range: -1
    java.lang.StringIndexOutOfBoundsException: String index out of range: -1
    String index out of range: -1
   </ows:ExceptionText>
  </ows:Exception>

There are no errors in the console.

What could be causing this? The ajax call is below:

var layerWFS = new ol.layer.Vector({
source: new ol.source.Vector({
    loader: function (extent) {
        $.ajax('http://servername:8080/geoserver/wfs', {
            type: 'GET',
            data: {
                service: 'WFS',
                version: '1.1.0',
                request: 'GetFeature',
                typename: 'lines',
                srsname: 'EPSG:2263',
            }
        }).done(function (response) {
            layerWFS
            .getSource()
            .addFeatures(new ol.format.WFS()
            .readFeatures(response));
        });
    },

    strategy: ol.loadingstrategy.all,
    projection: 'EPSG:2263',
    crossOrigin: 'anonymous',
}),

here is a snippet of the GeoServer log file:

2020-03-30 12:10:50,790 INFO [geoserver.wfs] -
Request: getServiceInfo
2020-03-30 12:10:50,806 INFO [geoserver.wfs] -
Request: getFeature
    service = WFS
    version = 1.1.0
    baseUrl = http://servername:8080/geoserver/
    query[0]:
        srsName = EPSG:2263
        typeName[0] = {http://boundlessgeo.com/}lines
    outputFormat = text/xml; subtype=gml/3.1.1
    resultType = results
2020-03-30 12:10:50,853 ERROR [geoserver.ows] -
java.lang.RuntimeException: java.lang.StringIndexOutOfBoundsException: String index out of range: -1
    at org.geotools.xml.Encoder.encode(Encoder.java:738)
    at org.geotools.xml.Encoder.encode(Encoder.java:609)
    at org.geoserver.wfs.xml.GML3OutputFormat.encode(GML3OutputFormat.java:302)
    at org.geoserver.wfs.xml.GML3OutputFormat.write(GML3OutputFormat.java:279)
    at org.geoserver.wfs.WFSGetFeatureOutputFormat.write(WFSGetFeatureOutputFormat.java:196)
    at org.geoserver.ows.Dispatcher.response(Dispatcher.java:996)
    at org.geoserver.ows.Dispatcher.handleRequestInternal(Dispatcher.java:279)
    at org.springframework.web.servlet.mvc.AbstractController.handleRequest(AbstractController.java:153)
    at org.springframework.web.servlet.mvc.SimpleControllerHandlerAdapter.handle(SimpleControllerHandlerAdapter.java:48)

here is what I found later in log under 'caused by':

Caused by: java.lang.StringIndexOutOfBoundsException: String index out of range: -1
    at java.lang.AbstractStringBuilder.setLength(Unknown Source)
    at java.lang.StringBuffer.setLength(Unknown Source)
    at org.geotools.gml2.simple.GMLWriter.coordinates(GMLWriter.java:302)
    at org.geotools.gml2.simple.GMLWriter.positions(GMLWriter.java:291)
    at org.geotools.gml2.simple.GMLWriter.posList(GMLWriter.java:348)
    at org.geotools.gml3.simple.LineStringEncoder.encode(LineStringEncoder.java:52)
    at org.geotools.gml3.simple.LineStringEncoder.encode(LineStringEncoder.java:34)
    at org.geotools.gml2.simple.FeatureCollectionEncoderDelegate.encodeValue(FeatureCollectionEncoderDelegate.java:196)
    at org.geotools.gml2.simple.FeatureCollectionEncoderDelegate.encodeFeature(FeatureCollectionEncoderDelegate.java:176)
    at org.geotools.gml2.simple.FeatureCollectionEncoderDelegate.encode(FeatureCollectionEncoderDelegate.java:130)
    at org.geotools.xml.Encoder.encode(Encoder.java:735)


texts:
So I use this query

            SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_dijkstra('
            SELECT gid AS id,
                     source::integer,
                     target::integer,
                     st_length(geom) AS cost
                    FROM leoforia',
            30, 40, false, false) ;

and as a result I get a table but when I try to open it in QGIS it says that table has no geometry. How to create geometry and generally visualize my result in OpenLayers?