texts: Format for viewing vector files offline and mobile

texts:
I wish to create a custom mobile application for Android and/or iOS for a school project that will allow a user to view (not edit, just view) specified ESRI shapefiles while offline.

Shapefiles at this link are an example of the source data I will need to use: http://www.fws.gov/southwest/es/AustinTexas/Maps_Data.html

What are good mechanisms for viewing Shapefile data offline on mobile devices? I would strongly prefer to use a format that is supported by free/open-source libs.

More specifically, is there FOSS library support for viewing the above example files directly on Android/iOS? Or do I need to convert the shapefiles to another format (I have heard of tiles, layers, geojson, database choices and so on) for viewing on mobile.

I would like to include these libraries in my own mobile app as I will be adding some features, i.e. not just view the files/data with another existing mobile app. I would also like to be able to include multiple shapefile sources as "layers" of a map in the app and ideally be able to interact with the layers as shapes. I.e. given a GPS point, I would like to find which polygons in included layers that point intersects. (Note: my question here is not really about how to do the point-in-polygon work, but rather to find a starting point for just viewing the Shapefiles in a mobile app that would enable later work)

Note: I am new to GIS so I am very new to the terminology and technology, but I have made a couple of android apps and some web apps so far.


texts: How to remove overlapping buffer in QGIS?

texts:
I wanna know how to do this in QGIS:

I've tried many ways, but still can't find the right solution.


texts:
I have been having trouble with MMQGIS geocoding that I am using behind proxy. Yes I had set the proxy settings from Settings>Options>Network correctly and was able to get connected to the internet but MMQGIS still didn't work and threw error upon running geocode. I could not find any useful help for this issue anywhere online, after weeks I realized the program works on another machine that had QGIS installed recently, this make me think there could be something about the version of QGIS, so I upgraded my QGIS to 2.14, but still the error persisted. Then it dawned to me to uninstall and reinstall MMQGIS plugin (yes even after uninstalling QGIS, the machine had preserved my older installation plugins). So I installed the recent version of MMQGIS and there we go, it worked. The issue seems to be the certain version of MMQGIS that I was using. I am copying the full error here in case this helps others, it wasted hours of my time.

Full error log with personal info masked:

  Traceback (most recent call last):   File
  "C:/Users/USERNAME/.qgis2/python/plugins\mmqgis\mmqgis_dialogs.py",
  line 1198, in run
      shapefilename, notfoundfile, fields, service, 1)   File "C:/Users/USERNAME/.qgis2/python/plugins\mmqgis\mmqgis_library.py",
  line 1890, in mmqgis_geocode_web_service
      x, y, addrtype, addrlocat, formatted_addr = mmqgis_geocode_address_google(address)   File
  "C:/Users/USERNAME/.qgis2/python/plugins\mmqgis\mmqgis_library.py",
  line 175, in mmqgis_geocode_address_google
      xml = urllib.urlopen(url).read()   File "~QGISsrc\apps\Python27\lib\urllib.py", line 87, in urlopen
      return opener.open(url)   File "~\QGISsrc\apps\Python27\lib\urllib.py", line 208, in open
      return getattr(self, name)(url)   File "~\QGISsrc\apps\Python27\lib\urllib.py", line 345, in open_http
      h.endheaders(data)   File "~\QGISsrc\apps\Python27\lib\httplib.py", line 969, in endheaders
      self._send_output(message_body)   File "~\QGISsrc\apps\Python27\lib\httplib.py", line 829, in _send_output
      self.send(msg)   File "~\QGISsrc\apps\Python27\lib\httplib.py", line 791, in send
      self.connect()   File "~\QGISsrc\apps\Python27\lib\httplib.py", line 772, in connect
      self.timeout, self.source_address)   File "~\QGISsrc\apps\Python27\lib\socket.py", line 571, in
  create_connection
      raise err IOError: [Errno socket error] [Errno 10060] A connection attempt failed because the connected party did not properly respond
  after a period of time, or established connection failed because
  connected host has failed to respond


texts: Why does Mapserver (ms4w) not generate tmp images?

texts:
I am having trouble with my web-based application, where my tmp images are not being stored for later retrieval. The images are generated each time they have to be loaded, so zooming out and then back in to the same location requires mapserv.exe to re-generate the tiles for that view.

The application lives here: http://rinnerweb.arts.ryerson.ca/moot2

A sample map file snippet is below:

    OUTPUTFORMAT                                        #-> BEGIN OUTPUT FORMATTING DEFINITION
        NAME "png"                                      # Name of format
        DRIVER AGG/PNG                                  # Image driver
        MIMETYPE "image/png"                            # Type of image
        IMAGEMODE RGBA                                  # Image mode
        EXTENSION "png"                                 # Image filetype
        TRANSPARENT ON                                  # Transparency
    END                                                 #-> END OUTPUT FORMATTING DEFINITION

    WEB                                                 #-> BEGIN WEB DEFINITION
        IMAGEPATH "/ms4w/tmp/ms_tmp/"
        IMAGEURL "/ms_tmp/"
    END                                                 #-> END WEB DEFINITION

The web page's source will show the OpenLayers init code, but I've pasted some sample layer init code below:

hii_1 = new OpenLayers.Layer.MapServer("Land Cover","../../cgi-bin/mapserv.exe",{
    map:'C:/ms4w/Apache/htdocs/hii/hii_landcover.map'},{
    isBaseLayer:false,
    transparent:true,
    format:"image/png",
    alpha:true,
    displayInLayerSwitcher:false
    });

Any ideas?


texts: Unable to use grass GIS algorithms in QGIS

texts:
I am unable to use GRASS GIS algorithms like v.distance in QGIS. Earlier I used to use the same from QGIS before opening standalone GRASS GIS GUI.
Is there anything I need to do to use GRASS algorithms in QGIS.


texts:
I am trying to get the area of land use in each pixel of the fishnet file. I used the following script:

arcpy.gp.TabulateArea_sa("fishnet_layer","FID","lulc.tif","Value","lulc_tab","1.1E-3")

Here fishnet_layer is the layer of fishnet shapefile which consists of 2000 pixels and lulc.tif is the land use image.

The above script gives the error:

ExecuteError: ERROR 999999: Error executing function.
("esriDataSourcesRaster.GdalRasterDataset") Failed to open raster dataset
ERROR 010419: Unable to build the attribute table of raster as required by the operation.
Failed to execute (TabulateArea).


texts: Incorrect raster extent in R but not QGIS?

texts:
I am attempting to read extents in R from a number of raster digital elevation models downloaded from here: http://www.tngis.org/elevation.htm

In doing so, I'm getting incorrect extents for a third (105) of the files in which the extents exactly match those of other files (alamo.dem=waterville.dem, armorel.dem=shooks_gap.dem etc.) :

library(raster)
extent(raster("alamo.dem"))==extent(raster("waterville.dem"))
[1] TRUE

In the above case, the extent of alamo.dem is correct, but the extent of waterville.dem is incorrect and should be more like (not exactly) (850440,861752,3963165,3977490). I know this because the waterville.dem correctly displays in QGIS even though its metadata extent is still wrong:

This must mean there is something in the file that correctly aligns the raster in QGIS other than the metadata and that R is not using this correct alignment information.

I'm wondering how to access this correct extent information to overwrite the incorrect extent information.

I've tried reading the binary data in R as well as opening the .dem in a word processor but I can't find anything that looks like the correct extent.

UPDATE:
as mdsumner suggested below,I tried opening waterville.dem in QGIS in its own map and got the expected, both the map extent AND the metadata extent were wrong.

When the map CRS is changed to UTM Zone 16N (EPSG 26916 or 26716), the map extent is correct for waterville. Thus, it seems that although the given crs for waterville is UTM zone 17N, its information is stored as if its in UTM zone 16N. Whoever created the dem incorrectly attributed it to UTM 17N when it should have been UTM 16N.

Although the stars package works as shown by mdsummer, the same can be done in raster:

library(raster)
r <- raster("waterville.dem")
r2 <- projectRaster(r,crs="+proj=utm +zone=16 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
extent(r2)
class      : Extent
xmin       : 849536.3
xmax       : 862083.7
ymin       : 3962859
ymax       : 3977785

This is the correct extent. Thus, the way to avoid this is to systematically check if each dem is in the correct CRS, and correct accordingly if not.


texts: Change background color of Layer in the TOC / QgsLayerTreeLayer with pyqgis

texts:
I would like to change the background-color (or font color) of some items in the TOC of QGIS:

Example-Mockup:

I tried several ways but nothing did the trick so far:

layer_list = QgsProject.instance().layerTreeRoot().children()
for lyr in layer_list:
    if isinstance(lyr, QgsLayerTreeGroup):
        lyr.setProperty("background-color","red")
        #lyr.setBackgroundColor(0, QtGui.QColor('green'))
        #lyr.setStyleSheet("background-color: green;")

    elif isinstance(lyr,QgsLayerTreeLayer):
        if lyr.name()=='Verteiler':
            lyr.setProperty("background-color","red")
            #lyr.setBackgroundColor(0, QtGui.QColor('green'))
            #lyr.setStyleSheet("background-color: green;")

I can't use a generic stylesheet like the ones used in the UI-themes, as the background-color should only be changed for some of my layers.

Is this possible?

edit1:

I found one way to change the style of layertree-items in general:

from PyQt4.QtGui import QApplication
QApplication.instance().setStyleSheet("QTreeView::item {color:green;}" );

Perhaps it would be possible to assign custom properties to some items of the TOC to be able to style only the wanted items.


texts: How to count the number of shapefiles that touch each polygon?

texts: I am trying to assign a number value field to different habitat shapefile polygons that is the number of other habitat type shapefiles that overlap or share a boundary. For example, I have 5 different habitat type shapefiles.  I want each habitat polygon within each habitat shapefile to have a value indicating the number of different habitats that it touches.  So if a particular wetland polygon touches a forest polygon and a beach polygon, that wetland polygon would receive a score of 2.  Thanks!


texts: Highlighting different features in layer with different colors in PyQGIS

texts:
I would like to highlight two features in a layer with two colors at the same time.
This is my highlight function:
def highlight_selection(self, color, results):
    """
    Highlight the selection on the canvas.  This updates all selected objects based on the result set.
    :param results: A dict-of-list of layer-features.
    """
    #self.clear_selection()
    for layer, features in results:#.items():
        print(layer)
        print(features)
        band = self.selectionbands[layer]
        print(band)
        band.setColor(color)
        band.setIconSize(25)
        band.setWidth(5)
        band.setBrushStyle(Qt.NoBrush)
        band.reset(layer.geometryType())
        band.setZValue(self.currentfeatureband.zValue() - 1)
        for feature in features:
            band.addGeometry(feature.geometry(), layer)
    #self.canvas.update()

when I use code below to highlight different features,
highlight_selection(mapwin, QColor(0, 255, 0), getFeatures(city1))
highlight_selection(mapwin, QColor(255, 0, 0), getFeatures(city2))

only city2 can be highlighted.
So, I'm confused.


texts: Showing legend with categories using a hexadecimal color scheme

texts:
I have a polyline vector layer with in one column the route number (of the bus) and in another the hexadecimal color (eg 'FF88AA'). Preview:

I know how to style the file but only by using the single symbology, found here:
How to style a vector layer in QGIS using hexadecimal color code stored in attribute?
Is there a way by using 'categories' or something as I need a legend with the route number and the matching color. Is this possible in QGIS?


texts:
Is there any simple way of excluding layer attributes from a QGIS tool output?
If I right-click --> Export, I get all layer attributes in the export, including those that are "Hidden" in the Attribute Form or Attribute table.
I can delete them from the output afterwards, but it can take some time (eg address data with 2.6 million records and only 3 or 4 of 85 data fields required).
It's the same with Dissolve results - the fields not selected to dissolve by are still present, but are worse than useless, because they have some random value in (presumably the from the first record encountered).
Ideally, it would be simple enough to give to those of my users who aren't up to writing SQL queries.
It was easy in ArcMap - get a list of field tickboxes with select all and unselect all buttons below, and anything unticked was omitted.
Solved:
I've been exporting by right-click --> Export --> Save features as, when I should be using Layer --> Save As


texts: Automatically update feature class based on another feature class in same geodatabase

texts:
I work at headquarters and there are 3 regions. Each region is responsible for their data and headquarters is responsible for merging the data and creating a national coverage. All regions and headquarters have the same fields in their feature classes.

I am using ArcGIS 10.1. The data is currently stored as follows: there’s one SDE geodatabase, each region has a Feature Dataset with their data in it. Headquarters also has a Feature Dataset with national data – which is just basically a merge of the regional datasets. What I want is an automatic or even semi-automatic way of updating a national feature class when the related regional feature class has been updated.

I originally wanted to have a geodatabase for each region and then create a parent child relationship but the IT people will not allow it.  We are allowed only one geodatabase on the network for this project.

I made a test file geodatabase locally and I have tried to create a Relationship Class with one region and headquarters. I put the region as the origin and national data as the destination. I made it a 1 to 1 and composite relationship. It was created successfully, but then when I edited the regional data (I added a new polygon, moved a polygon, deleted a polygon and changed the name of a polygon) the changes were not reflected in the related national feature class.

I want to know if the Relationship Class the best route to take for this project (assuming I can get it to work) or should I get my hands dirty with Python and use an update cursor? Or is there some other solution that I haven’t considered?

edit:
The database is a SQL Server Enterprise Database and Headquarters has all privileges, but the regions are only able to edit their own.  There are four feature classes per region, three of which will be edited. They are not currently versioned. Right now the process of updating is tedious – they send us the file and we update the national. They do all have GUID columns. Everything is running through shared services on Citrix, so there is a login associated with each user.


texts: I've received wind speed data in grid text format and I have no idea how to open it in QGIS. Please go easy on the computer jargon, I'm no computer scientist.


texts: Returning image for each iteration using Google Earth Engine?

texts:
I'm trying to map the maximal amount of days without rain using the TRMM dataset in Google Earth Engine. I'm doing this by iterating over the collection and if no rain has fallen, one gets added to the cell (code below). When rain has fallen the value gets multiplied by 0 so the "counter" is reset. Then I would like to store each result of each iteration in a Image collection and then select the maximum value to get the longest consecutive dry period.

But that is the theory. When I put this in a script I get an error while adding the image of one iteration to a list.

Does anyone know why this is and how this can be solved?

Code:

var list = [];
function drylength(current, previous){
  var mask = current.remap([0], [1], 0,"precipitation").rename('precipitation');
  var sum = mask.add(previous).multiply(mask);

  list.push(sum);
  return sum;
}

var dataset = TRMM
.filterDate("2015-01-01","2015-02-01")
.sort('system:time_start:');

var totalPrecipitation = dataset.iterate(drylength, dataset.max()
.remap([0], [0], 0, "precipitation")
.rename('precipitation'));

print(list);
print(totalPrecipitation);
Map.addLayer(ee.Image(totalPrecipitation), imageVisParam);

Furthermore it appears only 3 items are stored in the list which makes me assume the iteration is more complex than a literal iteration where all images are calculated one by one?

Here is an image of the error:

Errors written if image is not visible or for search engines:

Failed to decode JSON.
Error: Field 'value' of object '{"type":"ArgumentRef","value":null}' is missing or null.
Object: {"type":"ArgumentRef","value":null}.

and:

Unknown variable references: [_MAPPING_VAR_0_0, _MAPPING_VAR_0_1].


texts: Can I have two templates (content.ftl) on Geoserver?

texts:
I am using GeoServer for creating two web map applications. Under the GeoServer data directory, I have changed the content.ftl files for the use of my first map application. Now I am building the second map app, and would like to apply a different template to it.

Is it possible to define different templates for the two map applications?


texts: Comparing two shapefiles of US States using ArcMap

texts:
I would like to use ArcMap to create an overlap of two shapefiles; i.e. Alaska vs Texas state. link (source: alaska.org/how-big-is-alaska )

Is it possible to move the Alaska shapefile next and/or to overlap Texas state?


texts:
I'm using QGIS and I want to do a classification of an LandSat8 image.

I already uploaded the image and did a pansharping on it. I clipped the image to my area of interest. Then i set the bend set to get the different band combinations.

The problem is, when I choose 4-3-2 qgis projects 5-4-3 ... I sorted the single bands when setting the band composite to the right order.

So there must be a failure which i don't get.

I am working with the Semi Automatic Classification plugin.


texts: Variable substitution not working if multiple layers are provided in WMS call

texts:
I use cql expressions to make my color map entries dynamic as described here: http://docs.geoserver.org/2.8.x/en/user/styling/sld-reference/rastersymbolizer.html#cql-expressions

It works basically, but as soon as I add an additional layer in my WMS call (&layers=layer1,layer_raster) the parameters I pass in the URL (&env=low:10;medium:100;...) are just ignored and the default values I defined with the cql expression are applied.

Does anyone know about this issue? I am using GeoServer 2.8.0.


texts: Adding information from a polygon layer to my attribute table from point layer

texts: I have a point shapefile layer of animal observations and a polygon layer with 1 km x 1 km polygons of road density values. I'd like to add a column to the attribute table from the road density polygon data to my point layer. I'd like to grab the road density values that the observation points are located in and add them to my observation points attribute table. Do you have an idea of which tool I should be using for this?


texts: What is the correct tileInfo for MapBox basemaps in the ArcGIS JavaScript API (3.16)?

texts:
I'm trying to use a MapBox basemap in a JS API (3.16) app and they are using tiles that are 512px instead of the typical 256px. This causes the map to be one zoom level off and the labels and such appear 'shrunk':

Mapbox Satellite Streets | Mapbox
versus
Edit fiddle - JSFiddle

I think I might be able to set the tileInfo on the webTiledLayer to adjust for this, but after some playing around, I have not found success. Leaflet has a zoomOffset property in their tiledLayer that allows for this change, but no such thing exists in the ArcGIS JS API. I've also tried to modify the LODs, but I couldn't work either. What is the correct tileInfo for MapBox basemaps or how do I get the zoom levels adjusted so that the map renders correctly?


texts: Exporting county level data in *.shp file to text (csv) file in ArcMap?

texts:
I have a shape file created as per the answer to my previous question, which includes county level means for each county in the US.

How can I now export this data in a .csv file, e.g. with just the state, county, and county_mean fields?


texts: How can I stop QGIS from drawing a line when moving a label?

texts:
I have labels of a shapefile which I created through Easy Custom Labeling (labeling toolbar). I want to move the labels around manually to arrange them better but QGIS draws a line wherever I move them to. How can I disable this?

'Use custom dash pattern' was off anyway. I tried changing the line width to 0, even the pen to 'no pen' and also changed the value in the 'LblShowCO' column back to Null values. Still I have this annoying line :(


texts: Interface customization in QGIS -- I am unable to switch off the layer styling options of the layers

texts: I tried to switch off the layer styling using interface customization option but it is still visible in the layers tab. How to we restrict the styling as the project file will be used by many users .


texts: Trying to understand / troubleshoot intersection behavior of PostGIS

texts:
I'm trying to calculate the percent of intersection between two polygons using PostGIS.

I was able to do so successfully for 90% of my polygons, but about 1500 have very unexpected values.  I extracted one of these cases to examine the behavior.

The orange polygon resides w/in the purple polygons so the proportion should be near 100%.

However, after running:

UPDATE orange
SET pct_undevelopable = ST_Area(ST_Intersection(orange.geom, purple.geom)) / ST_Area(orange.geom)
FROM purple WHERE ST_Intersects(orange.geom, purple.geom);

I get the value for percent of area of intersection as: 0.000828452095134945

So I inspected the intersection part of the query and found that I got 3 rows as geometry outputs.

I exported all three as individual shapefiles: one simply doesn't work... I tried to open it in QGIS and ArcGIS and I also tried in CSV format for QGIS, no luck.

And two polygons.

Polygon 1 is what I imagine PostGIS is calculating percent of intersection from:

Polygon 2 is what I imagine I would like PostGIS to calculate it from:

I would love it if anyone can
1) explain what is going on (especially with the weird non functioning intersection geometry) and

2) how to fix it.


texts: FRAGSTATS Problem with "User Defined Tiles"

texts:
I have had recurring problems using FRAGSTATS 4.2.1 when trying to calculate four land cover classes over a large area using census tracts as the user-defined tiles. I have made sure the two layers I am using (land cover and census tracts, both converted to raster format) are identical in ArcGIS (e.g. column/row, uncompressed size, format (geoTIFF), pixel type, and pixel depth) and still have no success.

After running FRAGSTATS, it says "Error: Data buffering failure (level 2)" and "Analysis failed for file: tile_****" which that tile number is the Object_ID for the northernmost census tract in the study area. However, the output gives all of the patch, class, and landscape statistics for that "failed tile" and I have confirmed this by adding the "AREA" output column from the FRAGSTATS output and comparing it to the total area of the census tract in ArcGIS. So, why is it failing for that specific tile and only giving me the output for that tile, and why is it not giving me the output for my entire study area?


texts: Handling Nodata values in sklearn (working w pyimpute)

texts:
I'm using sklearn for a classification project. Because I'm relatively new to doing geospatial tasks outside of a GUI, I'm using Matt Perry's pyimpute (link has all the code for that module) to handle the heavy lifting of moving geospatial data into and out of sklearn-readable structures. It's built on rasterio and geopandas.

My ROI has a lot of zones that are irrelevant to the classification I'm trying to perform and I would prefer to just leave them out. If I include them in the classification, it will spuriously inflate the global accuracy of my model. But none of the materials I've managed to find on the topic address setting nodata masks for machine learning.

Is there a (relatively) simple way to accomplish this?

(I haven't posted any of my own code because this is a pretty general question)


texts: Any "standard" or at least popular API for geospatial data?

texts:
was looking at http://geogratis.gc.ca/api/en/ and had a brief look all over the internets. Came away wondering whether there is any form of mainstream open API used for providing access to geospatial data and perhaps metadata?

There's WFS of course but organizations appear content to build their own APIs. I haven't had a deeper look at parents/similarities of methods provided; thought I'd ask here first and see if anyone can give provide a general insight into the current scene.

Edit. In an attempt to narrow it down: If I had a database with spatial data that I would like to make available via an API, for integration in custom client applications, I would be joining others as this scenario appears to be happening in parallell in several places today - my question is if there's a more or less "standard" set of methods(provided by an API as described) or would a client application need to speak four different languages in order to query four different datasources via their different APIs, typically, right now?


texts:
I'm using Geowebcache seed command (via web interface at http://localhost:8080/geoserver/gwc/rest/seed/workspace:layerName) and it looks like everything runs smoothly but actually no tiles are saved, at least where I'd expect them to be saved (/data_dir/gwc).

I have GeoServer installed on top of Tomcat 7 as WAR. I also checked the logs (enabling GEOTOOLS_DEVELOPER mode) and it doesn't give any errors. Actually if enabling VERBOSE logging mode I also get something like...

...
2014-07-07 09:38:08,766 DEBUG [wms.map] - Writing png image ...
2014-07-07 09:38:08,771 DEBUG [wms.map] - Writing png image ... done!
2014-07-07 09:38:08,772 DEBUG [wms.map] - Writing png image ...
2014-07-07 09:38:08,776 DEBUG [wms.map] - Writing png image ... done!
2014-07-07 09:38:08,777 INFO [seed.SeedTask] - GWC Seeder Thread-1 completed (re)seeding layer workspace:layerName after 448 tiles and 45 seconds.
...

...and it look like everything went fine. But in the end I have no tiles.

By the way, the tiles I'm trying to save to EPSG:900913 projection coming from EPSG:3045 PostGIS store (already done this operation successfully on another VM having GeoServer installed on top of Tomcat6 via apt-get instead than via WAR).

Any suggestion on how to inspect/solve the issue?

Update:

Following @EvilGenius advice I was able to find the missing tiles (in fact they quite strangely were served by the Geowebcache to the frontend) in /tmp/tomcat7-tomcat7-tmp/geowebcache/workspace_layerName.

Now the question becomes: how comes they are not copied/moved to the final location (data_dir/gwc folder)? again, how can I inspect that?


texts: Add Leaflet layer based on variable with the same name?

texts:
I have a geojson layer named mapLayers and another geojson layer named lmiCounties, with a dropdown list with an item with the id lmiCounties. I want the second layer added to the first, which has already been added to the map when a user clicks on this item, and I want to do this with a function since I will have several layers/list items. The click capture is fine, its just that layer[0] (lmiCounties) does not associate with the actual lmiCounties geoJson layer. Seems like I am missing something simple here.

var lmiCounties = L.geoJson('some data');

document.getElementById('chooser').addEventListener("click", function(e) {
  mapLayers.clearLayers();
  layer = [e.target.id];
  console.log(layer[0]);
  mapLayers.addLayer(layer[0]);
  e.preventDefault();
});


texts: How to generate road curvature from the Points to path tool in QGIS?

texts:
I have been working on a project and I have used the "Points to path" tool in the processing toolbox to connect multiple points up along a route on the road network.

However, when processing this function it connects my points together but it only works in straight lines between sets of points. What I want is to have the lines between points follow the roads (routes) that would be used to reach the points.

I am currently running QGIS 3.2.2, I have tried searching the online forums to find the solution to no avail this may be down to my wording of the search.


texts:
I want to add some fine-grained noise to a raster in an automated python process (i.e. I do not want to use Photoshop or GIMP).  I looked at PIL, but unless I'm mistaken, it doesn't seem to have any modules for adding noise.  So, I'm thinking of using NumPy to create an array of Gaussian values in the three bands (RGB) and then add that to the raster (imported as an array to NumPy).

However, I'm sure somebody on the forum has done this sort of thing before and I imagine there could be better ways of going about it.

Does anybody have a recipe for adding noise to a raster using python (possibly with GDAL/Numpy/PIL/or something else).

Many thanks

EDIT:
Following Aragon's suggestion, I had a look at the noise module, but I need to use Python 3.2 and this module is quite old (2009) and is for Python 2.x.  So I decided to have a go at porting to 3 using the '2to3.py' tool.  All goes well until I run setup.py, where I get the following error:

running install
running build
running build_py
creating build
creating build\lib.win32-3.2
creating build\lib.win32-3.2\noise
copying perlin.py -> build\lib.win32-3.2\noise
copying shader.py -> build\lib.win32-3.2\noise
copying shader_noise.py -> build\lib.win32-3.2\noise
copying __init__.py -> build\lib.win32-3.2\noise
running build_ext
building 'noise._simplex' extension
creating build\temp.win32-3.2
creating build\temp.win32-3.2\Release
c:\Program Files\Microsoft Visual Studio 9.0\VC\BIN\cl.exe /c /nologo /Ox /MD /W
3 /GS- /DNDEBUG -Ic:\python32\include -Ic:\python32\PC /Tc_simplex.c /Fobuild\te
mp.win32-3.2\Release\_simplex.obj
_simplex.c
_simplex.c(224) : warning C4013: 'Py_InitModule3' undefined; assuming extern ret
urning int
_simplex.c(225) : warning C4047: '=' : 'PyObject *' differs in levels of indirec
tion from 'int'
c:\Program Files\Microsoft Visual Studio 9.0\VC\BIN\link.exe /DLL /nologo /INCRE
MENTAL:NO /LIBPATH:c:\python32\libs /LIBPATH:c:\python32\PCbuild /EXPORT:PyInit_
_simplex build\temp.win32-3.2\Release\_simplex.obj /OUT:build\lib.win32-3.2\nois
e\_simplex.pyd /IMPLIB:build\temp.win32-3.2\Release\_simplex.lib /MANIFESTFILE:b
uild\temp.win32-3.2\Release\_simplex.pyd.manifest
LINK : error LNK2001: unresolved external symbol PyInit__simplex
build\temp.win32-3.2\Release\_simplex.lib : fatal error LNK1120: 1 unresolved ex
ternals
error: command '"c:\Program Files\Microsoft Visual Studio 9.0\VC\BIN\link.exe"'
failed with exit status 1120

Any suggestions (either for an alternative library to try or how to port this one from 2 to 3)?  Should I take this to the main Stack Exchange site as we're drifting away from GIS now

Thanks


texts: Building customized data set from American Community Survey data?

texts:
I am trying to create a database in SQL Server with a customized dataset containing census information from the American Community Survey content available online. I want to show a granularity in my dataset and would like to download blocks & census tract level data of USA for 5 or 3-year timeframe ranging from [2009-2013]. I tried several methods to download both blocks & census tract level information. I could not get my hands on a data that would let me create tables and import values to them from the downloaded files.

I downloaded the block level data present in the PostgreSQL database section mentioned on the following link

http://www.justicemap.org/data.php

and got one 15GB file in the FILE format. Then, I tried to download a block group level data for USA from the link mentioned below

https://www.census.gov/geo/maps-data/data/tiger-data.html

and that gave me a collection of gdb related files. After which, I followed a couple of video and pdf tutorials to download content similar to this and got a lot of data with tons of label names with no reference. Could someone please help me obtain a good dataset in a usable format that contains block or census tract level information.


texts: How to create a histogram legend in QGIS print composer?

texts:
I've found this blog post while looking for a way to dynamically build a graph (histogram) in print composer.

The feature I want to replicate from the blog is drawing static rectangles in the qgis print composer and then size them based on a field value in this case as bars on a histogram.

Thus creating a legend/ histogram

I can follow this tutorial but I was just curious if a) there was an easier way outside of building the graph i need in matplotlib and importing an image or b) this functionality was going to be built in to qgis in future versions.


texts: How to export multiple MXDs to PNG without using VBA?

texts:
Okay, so I can feel that I am so close to solving this issue with the help of Jason (below) but I have modified to use as a script.  The GetParameterAsText returns the following error when trying to execute:

  : [Error 3] The system cannot find the
  path specified:
  u"'C:\Users\bob\Documents\folder\file.mxd';'C:\Users\folder\Documents\Folder\file.mxd';'C:\Users\folder\Documents\folder\file.mxd'\*.*"

Anyone know what is wrong with my code to keep getting this error?

    import arcpy, os, string

    #Read input/output folder path from script tool
    folderPath = arcpy.GetParameterAsText(0)
    outPNGfolder = arcpy.GetParameterAsText(1)

    for filename in os.listdir(folderPath):
        fullpath = os.path.join(folderPath, filename)
        if os.path.isfile(fullpath):
            basename, extension = os.path.splitext(fullpath)
            if extension.lower() == ".mxd":
                mxd = arcpy.mapping.MapDocument(fullpath)
                arcpy.mapping.ExportToPNG(mxd, basename + '.png')


texts:
I added the GeoBrain WCS server to the server list:
http://geobrain.laits.gmu.edu/cgi-bin/gbwcs-dem?version=1.0.0&

I connect to the server and get a list of layers. I select the first one listed - SRTM 30MUSA and select ADD. My map window is displaying a region in NM, the project CRS is EPSG:4326. The response is

Note that the bbox is nonsense.

Any idea what I'm doing wrong?

Edit:
This question is closely related, asking how to use a clip region to set the bbox in a WCS request, and the example is also GeoBrain.


texts: Creating plots and then extracting data from a TIFF file

texts:
So I have a giant .tif file with mean tree height data over the Everglades. I have a number of 10x10 field plots in the area.
For each plot, I have the center point and want to create a circle plot with a 5-meter radius and then clip and extract the data for the smaller plot from the larger TIFF file.

How do I go about doing this in QGIS? I used to know how to do it in ArcGIS but I can't run that on my Mac.


texts: Get Clipped geometry by a given bbox in postgis

texts:
I tried to get the geometry inside a bbox from postgis table, and I follow this instruction:

And I use the following sql ( the ? will be replaced by the bbox):

SELECT name, ST_Intersection(geom,ST_MakeEnvelope(?, ?, ?, ?, 4326)) as dest_geom  FROM  table_name where  st_intersects(geom,ST_MakeEnvelope(?, ?, ?, ?, 4326)) = true

However it seems that the returned features are not as expected, for example some features are not returned.

What's going  on?

And furthermore, as shown in the sql, the ST_MakeEnvelope(?, ?, ?, ?, 4326) occurs twice. Is it possible to avoid this?


texts: Map Design: Pros and Cons of using ArcGIS Halos to improve readability of text in a map

texts:
I have to create a map which involves 2 overlapping Layers. The underlying layer is coloured, the upper layer is hatched. On top of this I need to annotate the regions.

In ArcGIS the dirtiest method is to simply attach a halo to the text. Then the text is always readable. But I find the halo covers up so much area of the map and in print looks pretty awful.

I was thinking of exporting the geometries into The Gimp or Illustrator....but then my GIS Workflows become very complicated and updates are much more difficult to do.

Does anyone have any good examples of annotating under similar conditions and can offer some good advice?


texts: Specific methodology to call a GRASS function like v.transform with Python

texts:
I cannot understand if and how I can call with its parameters a GRASS function like v.transform through a Python script from e.g. PyCharm etc.
Is there something comprehensible to understand? An example would be very enlightening


texts: Dynamically adding items into one python addin Combobox from another?

texts:
I have two combo boxes. I am trying to add items dynamically into the second combobox from the first combobox's onSelChange(self, selection) event.
due to some reason, this is not working.

Here is my sample code:

import arcpy
import pythonaddins

class ComboBoxClass1(object):
    """Implementation for TestCombo_addin.combobox (ComboBox)"""
    def __init__(self):
        self.items = ["cb1item1", "cb1item2"]
        self.editable = True
        self.enabled = True
        self.dropdownWidth = 'WWWWWW'
        self.width = 'WWWWWW'
        self.cb2=ComboBoxClass2()

    def onSelChange(self, selection):
        pythonaddins.MessageBox(selection,"Message",0)
        self.cb2.items.append(selection)
        self.cb2.refresh()
    def onEditChange(self, text):
        pass
    def onFocus(self, focused):
        pass
    def onEnter(self):
        pass
    def refresh(self):
        pass

class ComboBoxClass2(object):
    """Implementation for TestCombo_addin.combobox_1 (ComboBox)"""
    def __init__(self):
        self.items = ["cb2item1", "cb2item2"]
        self.editable = True
        self.enabled = True
        self.dropdownWidth = 'WWWWWW'
        self.width = 'WWWWWW'

    def onSelChange(self, selection):
        pass
    def onEditChange(self, text):
        pass
    def onFocus(self, focused):
        pass
    def onEnter(self):
        pass
    def refresh(self):
        pass


texts: React-Leaflet-Draw: accessing a polygon's array of coordinates on save

texts:
I've got a component which puts an editable polygon on the map. When the user hits the "save" button, I want to access an array of the polygon's new vertices, so that I can save them.  How do I do this?

My component:

<FeatureGroup>
    <EditControl
        position="topright"
        onEdited={e => console.log(e)}
        edit={{ remove: false }}
        draw={{
                marker: false,
                circle: false,
                rectangle: false,
                polygon: false,
                polyline: false
             }}
        />
        <Polygon positions={polygonCoords} />;
</FeatureGroup>

The couple of references I've got:

https://github.com/alex3165/react-leaflet-draw

https://leaflet.github.io/Leaflet.draw/docs/leaflet-draw-latest.html#l-draw-event-draw:editstop

I understand I have to implement some sort of function dealing with the onEdited hook and the event generated thereby, but does anyone have any idea how I can get the new vertex array from this event?


texts: pyqgis: How to get features of selected geometries

texts:
I wonder if anybody could provide an example how to get the attributes of selected Features?

I tried the folowing Code in the python Console - but I'm stuck at the point where I'd like to get the attributes:

qgis.utils.iface
canvas = qgis.utils.iface.mapCanvas()
cLayer = canvas.currentLayer()
selectList = []
if cLayer:
    count = cLayer.selectedFeatureCount()
    print count
    selectedList = layer.selectedFeaturesIds()
    for f in selectedList:
        # This is where I'm stuck
        # As I don't know how to get the Attributes of the features

I would be really glad if anybody could help.

Regards


texts: Manually adding spatial reference information for FEWS rainfall data in ArcMap?

texts:
I am trying to match rainfall data from USGS/FEWS to coordinate locations of survey households in Uganda, so I will need to be able to project the rainfall data. So far, I downloaded the dekadal RFE 2.0 files for East Africa (bulk by year), and for each dekad, there are five files: .bil, .blw, .clr, .hdr, and .stx. When I add the data (the .bil file) to ArcMap 10.2.2, I get the warning:     "The following data sources you added are missing spatial reference information. This data can be drawn in ArcMap, but cannot be projected."

Can anyone advise on how to manually add the spatial reference for the data?

For reference, I downloaded the data files from: http://earlywarning.usgs.gov/fews/downloads/index.php?regionID=af&productID=3&periodID=6

The data documentation page http://earlywarning.usgs.gov/fews/africa/web/readme.php?symbol=rf, offers the following:

Dekadal RFE data characteristics:

Source: NOAA-CPC
Time step: 10-day (dekadal)
Resolution: 8km
Projection: Albers equal area conic
File Format: Windisp image, byte (8 bit); Generic BIL, "Integer" (16 bit).

Africa Continental Details:

 Coordinates for corners:

 Lower left lat : -42.243 deg
 Lower left lon : -23.490 deg
 Upper left lat : 43.711 deg
 Upper left lon : -24.600 deg
 Lower right lat : -42.242 deg
 Lower right lon : 63.414 deg
 Upper right lat : 43.712 deg
 Upper right lon : 64.523 deg

 Image size : 1152 rows x 1152 cols
 Center lat,lon : 1.000000, 20.000000
 Pixel size h x w : 8.000000 km x 8.000000 km
 Origin of latitudes : 1.000000 deg
 Central meridian : 20.000000 deg
 First std parallel : -19.000000 deg
 Second std parallel : 21.000000 deg
 Projection = ALBERS Conical Equal-area projection uses the clarke 1866 spheroid

I am happy to download the files for all of Africa instead of East Africa in case this somehow makes it easier to matching the documentation to the coordinates.


texts: View multiple mxd files in same tab of ArcGIS Pro

texts:
I have downloaded a bunch of map sheets which all fit together into one big map. I want to merge them all such that I can export any subsection of the combined map. They are provided as .mxd files and when I open them in ArcGIS Pro using 'Import map' each one opens in a separate tab. They are all georeferenced so it should be straightforward to merge them, but I can't find the option to do this anywhere.
What is the easiest way to achieve this?


texts: Reprojecting coordinates without building geometry object in ArcPy?

texts:
Is it possible with arcpy to reproject single x,y coordinates without building a complete geometry object?

For instance, something like this:

input = [[-8081441,5685214], [-8081446,5685216], [-8081442,5685219], [-8081440,5685211], [-8081441,5685214]]
output = [arcpy.A_Reproj_Fn((coord[0], coord[1]), source_SR, dest_SR) for coord in input]

I have huge arrays of coordinates to reproject, but for performance issue I do not want to build a geometry object before reprojecting.


texts:
Do you know of any free GIS sources of the Israeli administrative divisions: districts (mechozot, מחוזות) and subdistricts (nafot, נפות) which rely on official government information?

I heard there is a some free US government source, but failed to Google it.

Update (2014-02-11):

The Israeli government has opened a national GIS website. The UI will make your eyes bleed, the interface is only in Hebrew and there is no way to download the source files, only "save map as image".

On the other hand, it's the only official source. The map can be zoomed and downloaded chunk-by-chunk and georeferenced using QGis or alike. Depending on the required accuracy, this process can take anywhere between 30 minutes to a few hours and produce good-enough results.

Permalink to the Administrative boundaries layer.


texts: Stored Query on nested elements from an INSPIRE WFS 2.0 download service

texts:
Is there a way to create a stored query from this WFS 2.0.0 service that filters exactly an nested element?

http://geoportal.ancpi.ro/arcgis/rest/services/AU/AU_Download/GeoDataServer/exts/InspireFeatureDownload/service?VERSION=2.0.0&SERVICE=WFS&REQUEST=GetFeature&STOREDQUERY_ID=GETauByNationalCode&code=11940

This is an extract of the GML:

<au:AdministrativeUnit gml:id="auAdmUnitS.1387">
   <au:nationalCode>11940</au:nationalCode>
   <au:inspireId>
     <base:Identifier>
      <base:localId>1.29.11940</base:localId>
      <base:namespace>RO.Ancpi.AU</base:namespace>
      <base:versionId xsi:nil="true"/>
     </base:Identifier>
    </au:inspireId>
</au:AdministrativeUnit>

I would ideally like to be able to filter the data based on these two properties

au:InspireId/base:Identifier/base:namespace
au:InspireId/base:Identifier/base:localId

To do this I created a stored query with the following definition:

<DescribeStoredQueriesResponse xsi:schemaLocation="http://www.opengis.net/wfs/2.0 http://geoportal.ancpi.ro/arcgis/rest/directories/arcgisforinspire/schema/ogc/wfs/2.0/wfs.xsd">
    <wfs:StoredQueryDescription id="GetSQ2">
        <Title>GetSQ2</Title>
        <wfs:QueryExpressionText language="urn:ogc:def:queryLanguage:OGC-WFS::WFS_QueryExpression" returnFeatureTypes="gml:AbstractFeature">
            <wfs:Query typeNames="au:AdministrativeUnit">
                <fes:Filter>
                    <fes:And>
                        <fes:PropertyIsEqualTo>
                            <fes:ValueReference>au:inspireId/base:Identifier/base:localId</fes:ValueReference>
                            <fes:Literal>1.29.11940</fes:Literal>
                        </fes:PropertyIsEqualTo>
                        <fes:PropertyIsEqualTo>
                            <fes:ValueReference>au:inspireId/base:Identifier/base:namespace</fes:ValueReference>
                            <fes:Literal>RO.Ancpi.AU</fes:Literal>
                        </fes:PropertyIsEqualTo>
                    </fes:And>
                </fes:Filter>
            </wfs:Query>
        </wfs:QueryExpressionText>
    </wfs:StoredQueryDescription>
</DescribeStoredQueriesResponse>

When I create a GetFeature Request using this Stored Query:
http://geoportal.ancpi.ro/arcgis/rest/services/AU/AU_Download/GeoDataServer/exts/InspireFeatureDownload/service?VERSION=2.0.0&SERVICE=WFS&REQUEST=GetFeature&STOREDQUERY_ID=GetSQ2

I get this error message:
Operation request contains an invalid parameter value. The PropertyName "au:inspireId/base:Identifier/base:localId" is in step "basex:Identifier" not compliant with the Application schema.


texts:
What datum (reference ellipsoid) does Nasa Worldwind use?

I can't seem to find where Nasa Worldwind tells explicitly what datum they use. I need some sort of valid link.


texts: Remove artifacts from OSM roads data in QGIS

texts: I downloaded OSM data for Amsterdam from the bbike website. Opening the roads layer in QGIS I get these strange artifacts where road sections cross the whole map in a straight line. What would be the best way to correct/remove them? Has anyone else had this effect with OSM data?


texts:
What is a *.pkinfo file?

Google Search results talk about Package Info files and how they contain map documents and layers, so I don't see how this is different from a .mpk file.

I also have not found anything about how to create a *.pkinfo file.


texts: Adding to group layers in loop is slow in ArcGIS Engine?

texts:
I am using ArcGIS Engine 10.5.1.

I wanted to ask and see if anyone else has encountered this and or found a work around.  I have a loop in c# where I am taking a large set of layer files and adding them to the map.  I have tried this two ways but both ways seem to be equally slow (~10 minutes for ~180 layers).  The contents of the layers files are just shape objects.  Not a lot of data.  The first way I tried to achieve this is

IGroupLayer groupLayer = new GroupLayer();
foreach(SortedListLayers layer in MySortedListLayers)
{
    ILayerFile layerfile = new LayerFile();
    layerfile.Open(layer.location);
    groupLayer.Add(layerfile.Layer);
    layerfile.Close();
}

On the groupLayer.Add line there is more and more overhead  each time the call is made to complete adding the layerfile.Layer to the group even though the group is not yet on the map.  My alternative was to use IMapLayers and implement IEnumLayers on my SortedListLayers object (which only contains a sorted list of 'layers' and the location that it belongs to) and use the maplayers.InsertLayersInGroup call to try to diminish looping but that turned out even slower.


texts: postgis st_intersects returns false even though locating the address on google maps shows that it lies inside polygon boundary

texts:
I have an address in lat long format. When i pull up the address on Google Maps, the marker visually looks like it falls inside the map polygon boundary I have.

However, when I use the ST_INTERSECTS query to check the lat lon against the polygon boundary, I get a false. This is the query I am running

select st_intersects(st_geomfromtext('POINT(lon lat)', 3785), st_geomfromtext(polygon_coordinates));

I have checked the SRID of both the geoms to make sure they match and also tried ST_WITHIN and ST_CONTAINS. No luck.

Can someone throw some light on what I am missing here?


texts:
Recently I've reviewed ASTER GDEM v2 dataset and found many tiny (at least > 1ha) islets that's missing in other dataset such as WVS. I've change the altitude to 10m and they're still there.

Anyone/Organization ever attempted to documented/confirmed those? Or those are just defects?
I can't see them from Google Maps as the area isn't part of interest (no data).

Note:
- ASTER itself limit their scope (1 square degrees) to area with at least a significant island. So there are possibility of unknown islets everywhere else.
- The area that I reviewed is not area of polar ice extent.


texts: GDAL program produces only zeroes, regardless

texts:
I have been working on this problem for a while now and have asked two seperate questions regarding it (Counting events from a GRIB raster, ValueError: array larger than output file, or offset off edge Python GDAL), that have helped me tremendously. Still, the final solution eludes me. The code should iterate through a 24 images and compare each value to the conditions set in the loop, then write the number of instances the conditions were met into a GeoTIFF image. Despite being as sure as I can be, that this code should work as intended, I have still not been able to produce anything more than a GeoTIFF full of zeroes.

import gdal
import numpy
import osr
import os

del_mapa = str(input('Vpiši pot do datotek v obliki "C:\\Mapa1\\Mapa2\\...\\Zadnja mapa": '))
os.chdir(del_mapa)

pix_vel = 1000
#Defines pixel size
ime_rez = str(input("Poimenuj datoteko z rezultati: ") + ".tif")
#User can provide a name for the output picture
format = "GTiff"
driver = gdal.GetDriverByName(format)

vir_mere = gdal.Open('inca_20140101-0100.grb')
trans = vir_mere.GetGeoTransform()
geotrans0 = trans[0]
geotrans1 = trans[1]
geotrans3 = trans[3]
geotrans5 = trans[5]
#Extracts transformation data from one of the GRIBs
dst_ds = driver.Create(ime_rez, 401, 301, 1, gdal.GDT_CInt32)
dst_ds.SetGeoTransform([geotrans0, geotrans1, 0, geotrans3, 0, geotrans5])
#and pastes it directly into the newly created GeoTIFF
projInfo = vir_mere.GetProjection()
srs = osr.SpatialReference()
srs.ImportFromWkt(vir_mere.GetProjectionRef())
dst_ds.SetProjection(srs.ExportToWkt())

raster = numpy.zeros((301, 401), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)

def krog1(baza,novi):
    x = 0
    while x < 25:
        baza = gdal.Open(ime_rez)
        if x < 10:
            novi = gdal.Open('inca_20140101-%s00.grb') % ('0' + str(x))
            #This should open the next file in line
            if novi > 1 and novi < 12:
                baza = baza + 1
            else:
                baza = baza + 0

        elif x > 10:
            novi = gdal.Open('inca_20140101-%s00.grb') % (str(x))
            if novi > 1 and novi < 12:
                baza = baza + 1
            else:
                baza = baza + 0
        x += 1
    return dst_ds.GetRasterBand(1).WriteArray(baza)

dst_ds = None

Edit: The second half with additions and corrections as per the advice provided.

raster = numpy.zeros((301, 401), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)
banda = vir_mere.GetRasterBand(1)
bandtip = gdal.GetDataTypeName(banda.DataType)

def krog1(baza,novi,bandtip):
    vr_x = 0
    while x < 25:
        baza = gdal.Open(ime_rez)
        baza2 = baza.GetRasterBand(1)
        baza3 = baza2.ReadAsArray(0, 0, 301, 401).astype(numpy.CInt32)
        if x < 10:
            novi = gdal.Open('inca_20140101-%s00.grb') % ('0' + str(x))
            #This should open the next file in line
            re_raster = novi.GetRasterBand(1)
            band_raster = re_raster.ReadAsArray(0, 0, 301, 401).astype(numpy.bandtip)
            #Reads the raster as an array and then follows the condition checking
            if band_raster > numpy.bandtip(1) and band_raster < numpy.bandtip(12):
                baza = baza + 1
            else:
                baza = baza + 0
        elif x > 10:
            novi = gdal.Open('inca_20140101-%s00.grb') % (str(x))
            re_raster = novi.GetRasterBand(1)
            band_raster = re_raster.ReadAsArray(0, 0, 301, 401).astype(numpy.bandtip)
            if band_raster > numpy.bandtip(1) and band_raster < numpy.bandtip(12):
                baza = baza + 1
            else:
                baza = baza + 0
        vr_x += 1
    return dst_ds.GetRasterBand(1).WriteArray(baza)

dst_ds = None


texts: Merging 2 raster one with corn data one with wheat

texts:
This is my first time working with rasters outside of a class on remote sensing, Im working on learning the vocabulary and hope Im doing a ok job of it.

I would like to merge two rasters both have the same spatial extent one has number of harvested hectares for corn the other for wheat. Each raster cell has a value that gives the number of harvested hectares of corn or wheat. I would like the finished raster after the merge to show different color ramp for both corn and wheat. Each cells should show after the merge the value for harvested hectares area this could be corn or wheat.
I have used the raster calculator to merge the two raster but this gives me one value with a the two harvested hectares summed, I need have one raster with a column for corn harvested hectares and a column for wheat with harvested hectares.

I am using arcGIS 10.3


texts: Reprojecting layer in QGIS with python commands

texts:
I would like to reproject a layer from one coordinate system to another in qgis using python. The first answer to this post provides something similar to what I am trying to achieve, except that in my case:

I am using qgis 2.18
I would like to use actual layers instead of ``memory` layers and save the outputs to disk.

Here is what I wrote so far:

wkt = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'
in_crs = QgsCoordinateReferenceSystem(wkt)
wkt = 'PROJCS["WGS 84 / UTM zone 34N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",21],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32634"]]'
out_crs = QgsCoordinateReferenceSystem(wkt)

points = iface.activeLayer()
in_features = points.getFeatures()
xform = QgsCoordinateTransform(in_crs, out_crs)

out_layer = QgsVectorLayer('Point', 'points_out', 'memory')
out_features = []

#reproject
for feature in in_features:
    geometry = feature.geometry()
    geometry.transform(xform)
    feature.setGeometry(geometry)
    out_features.append(feature)

out_layer.dataProvider().addFeatures(out_features)

for f in points.getFeatures():
    print f.geometry().exportToWkt()

for f in out_layer.getFeatures():
    print f.geometry().exportToWkt()

#write out_layer
error = QgsVectorFileWriter.writeAsVectorFormat(out_layer, r'D:\11202750-002 - risicoanalyse en adaptatie experiment\Basis\point_new_crs.shp', 'CP1250', None, 'ESRI Shapefile')
iface.addVectorLayer(r'D:\11202750-002 - risicoanalyse en adaptatie experiment\Basis\point_new_crs.shp', 'point_new_crs', 'ogr')
print 'write file error = ' + str(error)

Which creates the following output:

Point (20.09595631708375407 39.81595508370384806)
Point (20.09814206032278605 39.81498734462326894)
Point (20.10044513203965622 39.81704382842367806)
Point (422623.74030901060905308 4407721.2567776320502162)
Point (422809.73621926747728139 4407611.96021181158721447)
Point (423009.15451480081537738 4407838.22327789757400751)
write file error = 0

Looking only at the console output it appears that the layer has indeed been reprojected. However, the layer that gets loaded by the last line is still in the original in_crs coordinate system:

And as it loads a warning is shown by the GUI:

How can I get the written file to be in the right CRS?

EDIT

I have also noticed that the "reprojected" point layer is placed at the right coordinates but in the wrong projection. That is the points are not in the same location. It seems to simply draw the points without taking into consideration that they are in a different coordinate system.

I have also tried replacing None in the last line by an actual QgsCoordinateReferenceSystem object as indicated by the API docs

error = QgsVectorFileWriter.writeAsVectorFormat(out_layer, r'D:\11202750-002 - risicoanalyse en adaptatie experiment\Basis\point_new_crs.shp', 'CP1250', out_crs, 'ESRI Shapefile')

but this throws a write_file_error = 6 and the same GUI warning as shown above. I'm still not sure why I can't set the coordinate system of my output file


texts: Leaflet: Style function does not work for point layer

texts:
I'm new to Leaflet and I'm not very proficient in Javascript. I've been trying to create a map visualizing a point layer with custom markers but for some reason the marker style as a function just doesn't work whereas marker style as a variable does. I've checked multiple resources and answers but I just can't seem to figure out why - it's probably something simple I'm missing. Here's the code that works:

var geojsonMarkerOptions = {
radius: 2.5,
fillColor: "#ff7800",
color: "#000",
weight: 1,
opacity: 1,
fillOpacity: 0.8
};

minesLayer= L.geoJson(
mines,
{ pointToLayer: function(feature, latlng){
    return L.circleMarker(latlng, geojsonMarkerOptions
    );
} }
).addTo( tnmap );

And here's the code that doesn't:

function minesStyle(feature) {
    return {
        radius: 2.5,
        fillColor: "#ff7800",
        weight: 1,
        opacity: 1,
        fillOpacity: 0.8
    }
}

minesLayer= L.geoJson(
mines,
{ pointToLayer: function(feature, latlng){
    return L.circleMarker(latlng, {
        style: minesStyle
        }
    );
} }
).addTo( tnmap );


texts: Shapefiles won't overlap due to projection problems

texts:
I have a two shapefiles. One of neighborhoods and one of points. The point file was originally a spreadsheet of X,Y coordinates. but have been exported to shapefile. They display but not over the neighborhood shapefile like they should. I've tried reprojecting,  project tool and define projection. It still will not overlap. Everything is displaying in US (Feet). What is wrong? Using arcmap 10.4

Here are the projection data in the Layer source:

Point file:
Projected Coordinate System:   NAD_1983_StatePlane_Pennsylvania_South_FIPS_3702_Feet
Projection: Lambert_Conformal_Conic
False_Easting:  1968500.00000000
False_Northing: 0.00000000
Central_Meridian:   -77.75000000
Standard_Parallel_1:    39.93333333
Standard_Parallel_2:    40.96666667
Latitude_Of_Origin: 39.33333333
Linear Unit:    Foot_US

Geographic Coordinate System:  GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree

Neighborhood file:
Projected Coordinate System:   NAD_1983_StatePlane_Pennsylvania_South_FIPS_3702_Feet
Projection: Lambert_Conformal_Conic
False_Easting:  1968500.00000000
False_Northing: 0.00000000
Central_Meridian:   -77.75000000
Standard_Parallel_1:    39.93333333
Standard_Parallel_2:    40.96666667
Latitude_Of_Origin: 39.33333333
Linear Unit:    Foot_US

Geographic Coordinate System:  GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree

Here's a sample of the Lat Long Coordinates:

Lat1        Long1
39.93128729   -75.14691545
39.93130836   -75.14711798
39.93143481   -75.14845446
39.93144498   -75.14853968
39.93117173   -75.14714819
39.93099668   -75.14705857
39.9310071    -75.14715449
39.93101477   -75.14725051
39.93119391   -75.14848865


texts:
I am aware that it is possible, in QGIS, to save your Vector Layer as a KML manually, however, I want to know if there is anyway I can do this programmatically from an external Python Script that I run from PyScripter.

I am aware ArcPy has that functionality, however, I am looking for a PyQGIS Solution.

Here is how I can do it with ArcPy using the Conversion Tool Set Layer to KML

...
Layer = Output_Layer
Output_File =  r"C:Users\path\KMZ.kmz"

# Process: Layer To KML
arcpy.LayerToKML_conversion(Layer, Output_File, "0", "false", "DEFAULT","1024", "96", "CLAMPED_TO_GROUND")
print "Layer to KMZ Conversion Complete"
print Output_Layer + " is now a .kmz"

or if there is a Python Library that will allow me to convert a shapefile to a KML/KMZ?


texts: How to add shapefile to dataframe from Python toolbox (ArcMap 10.4)

texts:
Does anyone know how I can add a shapefile to my current dataframe?

I have only been able to find examples that use a path for the shapefile.

    df = arcpy.mapping.ListDataFrames(mxd,"*")[0]

    PWshed = arcpy.RasterToPolygon_conversion("Wshed", "PWshed", "SIMPLIFY")

    arcpy.mapping.AddLayer(df, PWshed,"TOP")
    arcpy.RefreshActiveView()
    arcpy.RefreshTOC()


texts:
I am working in ArcGIS Desktop 10.5.1 (Basic level).

I have a point feature class and a polygon feature class. The point feature class has a field that contains either null or the ID of one of the polygon features. What's the easiest way to test whether each point is located inside the polygon it refers to and storing the result in an additional "inPolygon?" field on the point feature?

EDIT: I've got ~122k points and 33 polygons. Many of the polygons overlap. It's possible that a point could fall within all 33 polygons or within none of them, but each point references either none or exactly one polygon that I want to test for.


texts: Split a line feature using nearby points

texts:
I have a line shapefile and a point shapefile, and

the line has its own segments

and its attribute table,

the point shapefile has points near the line above, not on it

and the close look,

Can I use QGIS or other open source tools to split this line at locations which are nearest to the points？ And after this procedure I can get the original line with more segments(because being split). For example, find the location on the line, which is the nearest to the point, and split the line at this location(star).

I am using QGIS 2.2 under Windows 8.1 x64.


texts:
I have the following :

onEachFeature: function(feature,layer){
        console.log(this) // context is here(containing all needed f() & variables
        layer.on({
            click: function(e){

                }
         })
 }

I want to access the this context inside the click event click. I tried to bind(this) with layer.on.bind(this) but didn't work. So comes my question how can access to value stored in the higher context.


texts:
I have got a patch which supports Nvarchar and Nchar types when we need to publish sql server table to Geoserver. But I dont know, How to use this patch OR where i put this patch so that it can work with Geoserver ?

Can anybody explain, how i use this patch with Geoserver?


texts: How to cluster markers using the “Geojson-dashboard” in Leaflet JavaScript

texts:
I know how to use Leaflet.markercluster, but for some reason I can't use with geojson-dashboard framework template for Leaflet (http://fulcrumapp.github.io/geojson-dashboard/).

But I have a simple problem, I'd like to cluster the points but maintain all the setting, I've already tried so many options but until now I couldn't.


texts:
It has been a while since I have messed with the dataset options in the report designer. In the past I was able to set the dataset options to run reports on the current date and leave it like that within the report template.

The advantage to that is that anyone could run the report and it would generate results for that day.

However, I can't seem to remember the correct syntax. It was something like

Date = CURRENTDATE()

But that does not work and it will not verify.

Can someone please show me the correct way of doing this?


texts: Correct ASTER with SRTM or Topographic Map

texts:
I have a problem with my DEM again. I decided to use ASTER GDEM2 (for different reasons). ASTER GDEM2 have its problems, I know, and I dont want to discuss that again. I have to use it, because SRTM (USGS, CSI; Viewfinder) is not usable for my purpose.

So, there are 2 bad areas in my ASTER DEM:

These areas a nearly okay in CSI so I could fill these errors of the ASTER GDEM2 with it. But I guess, the result wont be as good as I want to (the mountain in fact is 500m from ground! SRTM say 350m). I also have a topographic map, a very old one, but the contours seems to be more accurate then SRTM and ASTER.

How can I correct the ASTER GDEM with SRTM or analog contours from my map?

Thanks!


texts: CartoDB Torque: show severity by size of marker?

texts:
I have a dataset containing the cost of damage caused by individual floods over the past 100 years.
I'd like to use TORQUE to display these events on a map over time, with their severity shown by the size of the marker. ie lesser cost has a smaller marker size, greater cost has larger size.

Is there a torque aggregation function which will allow me to do this?

Here is the style I have so far, but the events still show up as one size...

/** torque visualization */

Map {
-torque-frame-count:1024;
-torque-animation-duration:30;
-torque-time-attribute:"date";
-torque-aggregation-function:"count(cartodb_id)";
-torque-resolution:4;
-torque-data-aggregation:cumulative;
}

#cdd_cost_1{
marker-fill-opacity: 0.9;
marker-line-color: #FFF;
marker-line-width: 1.5;
marker-line-opacity: 1;
marker-placement: point;
marker-multi-policy: largest;
marker-type: ellipse;
marker-fill: #FF5C00;
marker-allow-overlap: true;
marker-clip: false;
}
#cdd_cost_1 [ value <= 956350000] {
marker-width: 25.0;
}
#cdd_cost_1 [ value <= 47368494] {
marker-width: 23.3;
}
#cdd_cost_1 [ value <= 21000000] {
marker-width: 21.7;
}
#cdd_cost_1 [ value <= 14523061] {
marker-width: 20.0;
}
#cdd_cost_1 [ value <= 8300000] {
marker-width: 18.3;
}
#cdd_cost_1 [ value <= 5738341] {
marker-width: 16.7;
}
#cdd_cost_1 [ value <= 4029790] {
marker-width: 15.0;
}
#cdd_cost_1 [ value <= 2828222] {
marker-width: 13.3;
}
#cdd_cost_1 [ value <= 1668590] {
marker-width: 11.7;
}
#cdd_cost_1 [ value <= 572824] {
marker-width: 10.0;
}


texts: Distinguishing between Layers and TableView using ArcPy?

texts:
I am working with ArcMap and am using ArcPy.

My task is to loop through mxds and find broken data sources.

My code is:

mxd = arcpy.mapping.MapDocument(fullPath)
brknMXD = arcpy.mapping.ListBrokenDataSources(mxd)
for brknItem in brknMXD:
    if brknItem.supports("workspacePath"):
        try:
            source = brknItem.workspacePath
            print  str(brknItem) + ": " +  source
        except:
            print "unexpected error"
            continue
    else:
        print "Layer does not support source"

Some mxds have TableViews and that is where the code breaks.

  Traceback (most recent call last):   File
  "P:\Scripts\Find_MXDs_Broken_Sources.py", line 5, in
      if brknItem.supports("workspacePath"): AttributeError: 'TableView' object has no attribute 'supports'

I checked within an mxd python interpreter to see what type of objects are return by arcpy.mapping.ListBrokenDataSources(mxd) and it indicated:

>>> for each in brknMXD:
...     print type(each)
...
<class 'arcpy._mapping.Layer'>
<class 'arcpy._mapping.Layer'>
<class 'arcpy._mapping.Layer'>
<class 'arcpy._mapping.TableView'>

Clearly I can just set a test to distinguish between these two types, but I don't know how exactly to write it.

Any ideas?


texts: Converting elevations into correct units with pyproj?

texts:
I am converting a point from wgs84 into a state plane coord system (epsg: 32051, which is defined in US survey feet) using pyproj, but the transform isn't changing the z values into feet:

point_4326 = (-81.23127933,37.70957554,668.105)
proj1 = pyproj.Proj(init='epsg:4326')
proj2 = pyproj.Proj(init='epsg:32051')
print pyproj.transform(proj1, proj2, *point_44326)

$ (589192.2346128735, 78765.53566868477, 668.105)

I have also tried specifying the exact proj4 string for the new coord system (which includes the USft scaling factor), but get the exact same result:

proj2 = pyproj.Proj('+proj=lcc +lat_1=37.48333333333333 +lat_2=38.88333333333333 +lat_0=37 +lon_0=-81 +x_0=609601.2192024384 +y_0=0 +ellps=clrk66 +datum=NAD27 +to_meter=0.3048006096012192 +no_defs')
print pyproj.transform(proj1, proj2, *point_4326)

$ (589192.2346128735, 78765.53566868477, 668.105)

Any ideas how I can get pyproj to convert the elevation values to the correct units, or if there are similar tools / libraries which can do this?


texts: Extracting vector information from a .dwg file into a QGIS environment?

texts:
I have a few CAD files in the .dwg format and I am trying to extract the data into a QGIS environment. The CAD files have a large number of different layers that represent different types of assets eg. water mains, underground electrical cables etc.
These layers are mostly polylines, but there are a few point features and polygons.

I believe the .dwg files have a geo-reference (there is an annotation that displays a Point Marker Number, Easting, Northing and AHD numbers)

Is there a way to accurately convert the polyline, point and polygon features from the CAD file into QGIS?


texts: Is there a better way to compare lists in arcpy?

texts:
I simply want to create a script which compares a list with the field names of a featureclass and if the fields do not exist in the featureclass add them.

My code works, but is very inefficient. How can I compare lists in a more efficient manner?

import sys,os
import arcpy
from arcpy import env

env.overwriteOutput = True

workspace = r"PATH_TO_WORKSPACE"
fc = r"PATH_TO_FEATURECLASS"

lst_fcfields = arcpy.ListFields(fc)
lst_myfields = ["field1","field2"]

for field in lst_fcfields:
        for myfield in lst_myfields:
                if field.name ==myfield:
                        print fc +" - same as - "
                else:
                        print field.name + " - not same as - " + myfield
                        arcpy.AddField_management(fc, myfield, "TEXT", "", "", 255)

Ideally, I would like the script to work interactively and ask the user if they would like to delete the fields which vary from the list, but as far as I know there are no dialogs in arcpy.


texts:
I'm trying to run a batch job for the Intersect(Analysis) tool in ArcMap 10.3, but when I do for 20 intersects, it brings back this error:

  ERROR 000210: Cannot create output.

However, when I do an individual intersect on the same shapefiles, it works perfectly fine.

I have tried restarting my session but get the same message, so is there I'm not aware of about this error?


texts: How to migrate PostGIS spatial data to MsSql spatial data?

texts:
am trying to migrate PostGIS geometry to SQL server Spatial data. which ETL i have to use. it's the first time i'll work with SQL server Spatial.
is there any Data ETL model that do this migration. thanks for help


texts:
I cannot export or add styles to AXF files, to style the AXF file in mobile device.

Using the ArcPad Check out Tool from Python, I can save an AXF file without styles because I reference the feature classes from the geodatabase, but I would like to add symbology style to this AXF file.

I tried exporting Stylesheet file .APS, I've added it to an Applet (copy and past XML text) but it didn't work. When I use the tool "ArcPad Check out" from ArcMap it exports with style inside every layer definition, I would like to know also if there's a manner to set Document (mxd) as Context in arcpy.


texts: QGIS 3 Display negative values bar diagrams

texts:
Suppose I have a dataset with an id (for each geographic entities) and 2 numeric fields NUM1 and NUM2. These numeric fields have positive and negative values.

When I use bar diagrams, only theses cases are displayed :

NUM1 is positive, NUM2 is positive
NUM1 is positive, NUM2 is negative
NUM1 is negative, NUM2 is positive

If NUM1 and NUM2 are negative, the bar diagrams are missing. Same story if we only use one numeric field. Is it a bug? How to display negative values in bar diagrams?


texts: Tilecache config for generation Mapnik tile

texts:
I have an XML stylesheet like below:

map:

<Map srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over" background-color="#b8dee6" >

A cover layer in:

<Layer name="cover" srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over">

and another like:

<Layer name="line" srs="+proj=utm +zone=39 +ellps=WGS84 +datum=WGS84 +units=m +no_defs">

And I can render picture with Mapnik via this script:

import mapnik
ll = (53.8889,26.481,54.0681,26.5971)
z=12
imgx = 50 * z
imgy = 50 * z
m = mapnik.Map(imgx,imgy)
mapfile = 't'
prj = mapnik.Projection("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over")
c0 = prj.forward(mapnik.Coord(ll[0],ll[1]))
c1 = prj.forward(mapnik.Coord(ll[2],ll[3]))
print c0,c1

mapnik.load_map(m,mapfile)
bbox = mapnik.Envelope(c0.x,c0.y,c1.x,c1.y)
m.zoom_to_box(bbox)
mapnik.render_to_file(m, 'test10.png')

But I can't generate tile with tilecache:

tilecache.cfg:

type=Mapnik
mapfile=/home/test.xml
srs=EPSG:900913
bbox=53.8944,26.4911,54.054,26.5820


texts:
When deriving NDVI out of a Sentinel-2 image - S2A_MSIL2A_20190502T074621_N0211_R135_T38SLF_20190502T102541 - I am getting a strange result when doing it in python vs in QGIS.

In QGIS (using raster calculator) the range of values is as expected, between -1 and +1, which I can confirm from the product metadata:

When doing the NDVI calculation in python usign the following code:

NIR = rasterio.open(S_files[3]).read() # band 8
Red = rasterio.open(S_files[2]).read() # band 4

NDVI_2 = (NIR - Red) / (NIR + Red)

I am getting the following min and max values:

NDVI_2.min()
Out[756]: 0.0

NDVI_2.max()
Out[757]: 21845.0

If you have a look the the histogram of NDVI_2

In addition, looking for the pixels that have value > 1:

print("Values bigger than 1 =", NDVI_2[NDVI_2> 1])

Values bigger than 1 = [ 4.7842736   4.72178932  4.71683968 ... 20.07305095 15.2136194
     17.2828496 ]

print("Their indices are ", numpy.nonzero(NDVI_2 > 1))

Their indices are  (array([0, 0, 0, ..., 0, 0, 0]), array([    0,     0,     0, ..., 10979, 10979, 10979]), array([ 2382,  2383,  2384, ..., 10444, 10445, 10446]))

Any idea what I am doing wrong?

--- EDIT ---

Testing the masked. vs non-masked statistics:

NON MASKED

red = rasterio.open(S_files[2]).read()
nir = rasterio.open(S_files[3]).read()
ndvi = (nir - red) / (nir + red)
ndvi.max()
Out[111]: 21845.0
ndvi.min()
Out[112]: 0.0

MASKED

red = rasterio.open(S_files[2])
nir = rasterio.open(S_files[3])
redc, out_transform = rasterio.mask.mask(red, features, crop=True, all_touched=True)
nirc, out_transform = rasterio.mask.mask(nir, features, crop=True, all_touched=True)
ndvi = (nirc - redc) / (nirc + redc)
ndvi.max()
Out[118]: 63.337864077669906

ndvi.min()
Out[119]: 0.0


texts: How to convert bearing, distance, xy coordinate measurements to point shapefile?

texts:
I have a natural resources dataset where individual seedlings in canopy gaps were recorded from a center "legacy tree" within each canopy gap.  Each sapling has the following measurements associated with it:

Legacy tree GPS xy coordinates (i.e. a center location in the canopy gap)
Bearing to sapling from the legacy tree (e.g. 218 degrees)
Distance in meters from the legacy tree

There are thousands of these individual tree locations, so an automated approach would be best.

How can I convert this dataset into a point shapefile/featureclass?  ArcGIS, Arcpy/Python and R solutions would be especially useful.


texts: How to select points within polygon constrained by field?

texts:
I've created a polygon from a raster dataset:

arcpy.RasterToPolygon_conversion(RasterInt, "C:/arc/poly_int.shp", "NO_SIMPLIFY", "Value")

The 'value' field has numbers 1-37, which map to the ID field in the output polygon shapefile feature class.

Now I would like to select points that lay within certain polygons (e.g. delineated by ID). Using the SelectLayerByLocation_management tool only selects points inside all the polygons. How to a define a criteria?

Alternative approach: How do I create multiple polygons from the main polygon, using the ID field?

I'm looking for solutions in Python -- not necessarily using arcpy but I do have it.


texts: Is it possible to use a local variable in a label expression?

texts:
Is it possible to use a local variable in a label expression/field calculator in QGIS 3? I have an extremely long CASE expression that generates a String for the labels of my objects, and I would like to process the string after I have generated it. To do this, however, I would need to save the string in a local variable and refer to the variable in the processing afterwards.

Here's a shortened version of my problem:

CASE when ("A" != 'null'  then 'A: ' ||  "A" || ', ' else '' end ||
CASE when ("B" != 'null' then 'B: ' ||  "B" || ', ' else '' end ||
CASE when ("C" != 'null' then 'C: ' ||  "C" || ', ' else '' end

This creates something like: "A: 3, B: 1, "

There's ", " always after the last variable and I would like to remove that from the end. The expression I'm using is so long that it exceeds the character limit of the field calculator. That's why I would like to save the string in a local variable "input", then use left(input, length(input) - 2). This would delete the last two ", " and make the string "A: 3, B: 1". It will have to be done like this because the table will be updated in the future.

Is this possible?


texts:
I'm trying to add a field to a created feature class, then create a key based on whether the AccountNo is equal to the imputed Feature1 through 6. For example, if Feature1 is equal to the AccountNo when it goes through the list of polygons, it will set the feature class in the created field to 1, if Feature2 is equal to the Account No, it will set it to 2, and so on. My code only has one update query now, but I'd like it to run through each "Feature". Is this doable using this technique?

import arcpy

FClass = arcpy.GetParameterAsText(0)
Field = arcpy.GetParameterAsText(1)
Feature = arcpy.GetParameterAsText(2)
Feature1 = arcpy.GetParameterAsText(3)
Feature2 = arcpy.GetParameterAsText(4)
Feature3 = arcpy.GetParameterAsText(5)
Feature4 = arcpy.GetParameterAsText(6)
Feature5 = arcpy.GetParameterAsText(7)
Feature6 = arcpy.GetParameterAsText(8)

relateFieldsList = ["StrNum", "Street", "StrSuf", "StrUnit", "AccountNo"]

for fc in arcpy.ListFeatureClasses("Comps"):
    arcpy.AddField_management(fc, "Comp_Num", "TEXT", field_length = 10)
    with arcpy.da.UpdateCursor(fc, "Comp_Num") as cursor:
        for row in cursor:
            if "AccountNo" == Feature1:
                row[0] = 1
            cursor.updateRow(row)


texts: Is there an equivalent of 'zonal statistics' (always used with a raster) for polygons?

texts:
I want to extract statistics from a polygon using other smaller polygons as the mask. The smaller polygons have unique id numbers that I want to assign the calculated statistics too.
If I wanted to get statistics from a raster, I would use zonal statistics and use my small polygons to delineate the boundary of where I want my statistics from in the raster (as in this post Iterating through multiple polygons to produce a table with elevations using DEM?).
The zonal toolbox only seems to accept rasters as the input value (see image below)

Is there a way that I can extract the statistics from a polygon (instead of a raster) using other polygons as the 'mask' or boundary of where the data is taken from?

I tried an intersect, but that just cut out the shape of the smaller polygon, without retaining all the attributes from both layers (I lost my unique id).

I also tried a spatial join, which seems to have pulled the data together from both attribute tables, but I don't think it actually calculated statistics on any of the values (I think it just returned one value from within the polygon boundary).

My smaller polygons are overlapping, so I also may need to split them first and iterate through whatever steps are suggested to me in ModelBuilder (as is the case in this post:Iterating through multiple polygons to produce a table with elevations using DEM?).

This is a screenshot of the output I get from 'intersect'...note that my stationnrs are all the same...as opposed to acting like unique identifiers


texts: Removing black pixels of georeferenced image without affecting interior pixels using GDAL?

texts:
I am looking for an algorithm or function in GDAL that can remove black pixels out of image boundaries. Note that I've tried gdal_fillnodata (gdal.FillNoData) and the result was good in most cases; But in some images I can find a lot of pixels in middle of image which are black and so they will be removed like out of boundary pixels.

What can I do for making difference between interior and exterior pixels and then removing them?

I use gdal with python.


texts: Coordinate Reference System change problem in GeoTools

texts:
I try to convert some coordinates from ETRS89-LAEA (EPSG:3035) to WGS84 (EPSG:4326) with geotools. Here is the code:

Coordinate c = new Coordinate(3847988.7, 3161178.2);
MathTransform mt = CRS.findMathTransform(
   CRS.decode("EPSG:3035"),
   CRS.decode("EPSG:4326"));
c = JTS.transform(new GeometryFactory().createPoint(c), mt).getCoordinate();
System.out.println(c);

As a final result, I get:

(56.34206943661075, -8.998327460134755)

which is obviously wrong. What is wrong with my code?

Note:

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

returns:

PROJCS["ETRS89 / LAEA Europe",
  GEOGCS["ETRS89",
    DATUM["European Terrestrial Reference System 1989",
      SPHEROID["GRS 1980", 6378137.0, 298.257222101, AUTHORITY["EPSG","7019"]],
      TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
      AUTHORITY["EPSG","6258"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic latitude", NORTH],
    AXIS["Geodetic longitude", EAST],
    AUTHORITY["EPSG","4258"]],
  PROJECTION["Lambert_Azimuthal_Equal_Area", AUTHORITY["EPSG","9820"]],
  PARAMETER["latitude_of_center", 52.0],
  PARAMETER["longitude_of_center", 10.0],
  PARAMETER["false_easting", 4321000.0],
  PARAMETER["false_northing", 3210000.0],
  UNIT["m", 1.0],
  AXIS["Northing", NORTH],
  AXIS["Easting", EAST],
  AUTHORITY["EPSG","3035"]]

and

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

returns:

GEOGCS["WGS 84",
  DATUM["World Geodetic System 1984",
    SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
  PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
  UNIT["degree", 0.017453292519943295],
  AXIS["Geodetic latitude", NORTH],
  AXIS["Geodetic longitude", EAST],
  AUTHORITY["EPSG","4326"]]

which both seem correct...


texts:
I installed QGIS 2.16.0, and for a while everything was ok. When I uninstalled some plugins, appeared to me the window, minidump.

I installed qgis 2.16.2. the problem still exists. when I install and uninstall any plugin I'm getting message box, dump, and QGIS close down. Where is the problem?

I still have the same problem with the dump file. For all versions of QGIS (14; 16; 18) again the same problem. This is happening on three different computers.

What is the problem?

This dump file has a few megabytes (33MB) and quickly consumes the memory. After a while I have to manually look for these files to be deleted.


texts: Spring and OpenLayers conflicts on label rendering?

texts:
I need to use Jquery to handle other elements on the webpage. I tried to put labels on a vector layer with the following style:

  var vector_style = new OpenLayers.Style({
            'fillColor': 'red',
            'fillOpacity': .8,
            'strokeColor': '#aaee77',
            'strokeWidth': 3,
            'pointRadius': 8,
            'label': "${name}",
            'labelAlign': "cc",
            });

But the labels do not show.

I am just wondering if there is a conflict on this one  "${name}".
I tested it again with a simple page without Jquery. The labels do show.
Please let me how to fix this. Thanks.

Update:
I just figured out the conflict may come from Spring or JSTL. I use Firebug tool to check the javaScript source code of the style. It appears as this:

    var vector_style = new OpenLayers.Style({
            fillColor: 'red',
            fillOpacity: .8,
            strokeColor: '#aaee77',
            strokeWidth: 3,
            pointRadius: 8,
            label: "",
            labelAlign: "cc",
            });

I returned an variable from the server side:

    model.addAttribute("name", "test");

then all features show the label "test". It approved that ${name} is recognized as a spring variable.


texts: Lambert93 to WGS84, IGN algorithms, longitude problem

texts:
I'm looking for a Java-based solution to convert coordinates from Lambert93 to WGS84.

I found this, "lambert-java", which is an implementation of the algorithms (PDF) published by IGN (French government agency in charge of geographic information).

However, when testing it with data from North of France, the resulting longitude is off by something like 50km. Here's an example:

Source data (in Lambert93): X = 668832.5384, Y = 6950138.7285
Geofree conversion to WGS84: lat = 49.64961, lng = 2.56865 (this is correct)
Output of the above Java program: lat = 49.65262, lng = 1.90503 (this is incorrect)

As far as I can tell, the program is "faithful" in that it faithfully implements the algorithms from IGN; so the fault either lies in the algorithms themselves, or in some parameter or some constant that needs to be adjusted.

How can one sort this out?

(There are other Java-based solutions capable of doing all sorts of GIS manipulations, but since I only need to do this conversion (on lots of data), I'd rather use a very simple, very light package, and "lambert-java" fits this bill nicely; however of course it also needs to be correct...)

Edit 2014-03: this was solved by the original author five months ago and has been working wonderfully ever since.


texts: Getting mean values from multiple extracted layers

texts:
Bit of a novice in using ArcGIS

I've used "Extract by Mask" -> Batch to get multiple layers clipped to a catchment area I have drawn as polygons.

Is there a way to get the mean values for each of the extracted layers at once?

Without having to go through each one individually?

At the moment I'm going through the properties of each individual created layer and getting the mean value from there. The issue is that I have to do hundreds of these.


texts: How to find a route between two points if the two are not along the network?

texts:
I have browsed through a lot of posts and I haven't seen any to address this issue.

My research plans to solve the distance from a point to another point using the selected highway network (selecting out other roads system). But I have points that are not on the highway system. I am wondering if there is any tool that can pinpoint reachable driveways that capture all points of interest (thousands) within a certain distance and calculate the distance between each two driveways?

Alternatively, let's say if I put all drivable roads on a map. How do I order ArcGis to calculate routes of points to the nearest highway node and routing to the destination merely choosing the highway instead of other roads?


texts:
I know how to get records in a shapefile using PyShp, but how to sort records by an attribute (i.e., NAME or POPULATION)?
My code:
import shapefile

file = "c:/data/cities.shp"

with shapefile.Reader(file) as shape_file:
    feats = shape_file.shapeRecords()

for feat in feats:
    print(feat)


texts: What iOS (iPhone/iPad) app from the App Store can I use to collect data?

texts:
What iOS app from the iTunes App Store can I use to collect data (mainly points, perhaps photos) and later import them into my GIS?

I am mainly interested in an out-of-the-box solution, not develop one using any API. Bonus if it already connects to either my own WMS, WFS or Tile server.

Update: Free or paid doesnt matter. I just want the best option.

Update years later!: Thank you all for all your answers. At that time, I needed to record points/lines/polygons and support workflows with relationships (e.g. if I needed to add signs to a light pole, they would not go in the light pole features but in a related table called signs). I ended up going through this rabbit hole, that once solved, was the genesis of the AmigoCloud mobile client (Open Source).

Here is a video of how it works: https://www.youtube.com/watch?v=lj6OeTZgzM8


texts:
I have a added Graticule control in my OpenLayers map. Now I want to toggle this using a control button. This what I have written:

function init() {
      graticuleCtl = new OpenLayers.Control.Graticule({
                    numPoints: 2,
                    labelled: true,
                    displayInLayerSwitcher: true
      });

    map = new OpenLayers.Map('map', {controls: []});
    map.addControl(graticuleCtl);
.
.
.
      var GraticuleCtl = new OpenLayers.Control.Button({
            displayClass: "olGraticule",
            title: "Show Graticule",
            trigger: showGraticule});

        // Panel (toolbar)
        Panel = new OpenLayers.Control.Panel({
            defaultControl: DragPanCtrl
        });
        Panel.addControls([
            .
            .
            .
            .
            GraticuleCtl
        ]);
        map.addControl(Panel);
}

function showGraticule() {
   if (graticuleCtl.visible==true){
       graticuleCtl.deactivate();
       return true;
       }
   else if (graticuleCtl.visible==false){
       graticuleCtl.activate();
       return true;
       }
   else{}
}

But this works for once only, since the 'visible' property checks only for initial map load.


texts: translating and scaling viewport to fit geojson map in D3js

texts:
I am trying to load a map of Rio de Janeiro (Brazil) in GeoJson for with D3js

the code and the map is here:
https://gist.github.com/4157853

I can open the map normally with Qgis, but with d3js nothing appears in the SVG canvas.
My projection is epsg4326, the same as google maps. and the units are in lat/long.

What am I doing Wrong?


texts: Seeking geoprocessing tool for adding vertex at intersection without splitting

texts:
I have a polyline shapefile and a polygon shapefile. I need to add vertices to the polylines where they cross polygon boundaries. The polyline length must stay the same, so the Intersect tool doesn't work.

What tool should I use? I just want to add vertices and I will be performing this for dozens of datasets.

EDIT: Below are two attempts that don't work. Code & expected behavior identified.

Base-Case to Check Expected Behavior:

outpath = r'C:\myfolder'
shp_input_lines = r'C:\myfolder\original_lines.shp'
shp_input_polygons = r'C:\myfolder\original_polygons.shp'

# make sure there are no multipart gremlins hidden in here
clean_lines = outpath + '\\' + '_clean_lines.shp'
arcpy.MultipartToSinglepart_management(shp_input_lines, clean_lines)

# add "my_id" field
arcpy.AddField_management(clean_lines, "my_id", "SHORT")
counter = 1
with arcpy.da.UpdateCursor(clean_lines, ["my_id"]) as cursor:
    for row in cursor:
        row[0] = counter
        cursor.updateRow(row)
        counter += 1

output1 = outpath + '\\' + '_clean_lines_intersected.shp'
arcpy.Intersect_analysis([clean_lines, shp_input_polygons], output1)

output2 = outpath + '\\' + '_clean_lines_intersected_vertices.shp'
arcpy.FeatureVerticesToPoints_management(output1, output2, point_location="ALL")

CHECKS, this is basis for expected behavior:

number of lines in "clean_lines" = 3123
number of points in "output2" = 14173

Try 1, Failed:
Intersect then Dissolve

# using output1 from above Base-Case
output3 = outpath + '\\' + '_clean_lines_intersected_dissolved.shp'
arcpy.Dissolve_management(output1, output3, dissolve_field="my_id", multi_part="SINGLE_PART", unsplit_lines="DISSOLVE_LINES")
# also tried with unsplit_lines="UNSPLIT_LINES" made no difference

output4 = outpath + '\\' + '_clean_lines_intersected_dissolved_vertices.shp'
arcpy.FeatureVerticesToPoints_management(output3, output4, point_location="ALL")

CHECKS:

number of lines in "output3" = 3123 (good)
number of points in "output4" = 13846 (bad; expected 14173)

Try 2, Failed:
Integrate polyline and polygon files

# using clean_lines from above Base-Case
cleancopy = outpath + '\\' + '_clean_lines2.shp'
arcpy.CopyFeatures_management(clean_lines, cleancopy)
arcpy.Integrate_management([cleancopy, shp_input_polygons], "0.1 feet")

CHECKS:

number of lines in "cleancopy" = 3120 (bad; expected 3123)


texts: Question about Using r.mapcalculator in Python console

texts:
According to the following conversation:
Using r.mapcalculator in Python console?

I tried:

##NoData Background=name
##RASTERLAYER_BAND=raster
##nodatabackground=output raster

import processing

ext = RASTERLAYER_BAND.extent()
a = str(ext.xMinimum())
b = str(ext.xMaximum())
c = str(ext.yMinimum())
d = str(ext.yMaximum())
size = a + ',' + b + ',' + c + ',' + d

nodatabackground=processing.runalg('grass:r.mapcalculator', RASTERLAYER_BAND,None,None,None,None,None,'A/A*A',size,0.0,nodatabackground)

But it doesn't work to me. I get the following message: AttributeError: 'unicode' object has no attribute 'extent'

Any idea to solve it?


texts: QGIS 1.8.0-Lisboa DBmanager (OSx) can't connect to remote PostGIS 2.0.2 database (Debian)

texts:
The PostGIS database is remote but I have made changes to the Postgresql config files and can access it remotely via pgadmin.

When I open DBmanager in QGIS I see a note on the bottom that says "postgis:No module named psycopg2." I can't find anyway to connect to the database through the GUI.

I have been able to connect to the database in QGIS using the SPIT plugin, though when I try to import layers it returns the following error: "Problem inserting features from file /path/to/shape/file.shp Invalid table name."

I can get the data from QGIS to the server by ftping the files to the server and then importing them with the shp2pgsql command. I can then download the data back into PostGIS without a problem, so I know the database is working properly.

What I want though is to be able to import QGIS layers to the PostGIS database from within QGIS.

Please forgive my ignorance. I'm learning this stuff as I go. Any clues you can provide will be much appreciated.

I've looked at the answer here: psycopg2 error with QGIS 1.5. I might be missing something, but it didn't help me. That answer explains how to config Windows registry keys, but I'm running QGIS on OSX and the database server is running on Debian (stable), so I don't have registry keys.


texts: How to export a map with a transparent background in QGIS?

texts: I have gone into Composer and specified the page background to be transparent.  It still comes out white.  Is there a solution?


texts: MapBox geocoder results are not exact, how to enhance mapbox geocoding API?

texts:
We are trying to get the address related data (ex. site address, city, state, country, postal code, lat-long)
of a particular site from the MapBox geocoding API, through geocoding and reverse-geocoding and populate it into form fields. We are using the autocomplete address search field for the same. but it is observed that whenever the user tries to search few sites with the help of autocomplete address search field the geocoding API does not return exact address or even related address into its suggested address dropdown (especially for sites in Poland, UK).

It will be great if someone could suggest any solution to enhance the geocoding API functionality so that a user can search address of any site
can we also design our own dataset and use it into geocoding API to enhance the API functionality. so that all the required
site address will be readily available for the user whenever he/she uses the autocomplete address search field (geocoding API)