Friday 30 June 2017

qgis - Converting from QML to SLD on the command line


We have the problem of automating the conversion of MapInfo TAB styles to SLD. I've managed to achieve the first part of the process (TAB -> QML) using this tool:


https://github.com/NathanW2/MapInfo-to-QGIS-style-generator


I can load the resulting QML file into QGIS and export it as SLD manually. However, would I really need is a way of doing this via the command line.


Is this possible with QGIS, or does anyone know another tools that might help here?





Installing PostGIS on older homebrew versions of PostgreSQL (e.g. postgresql@9.6)?



As there are two ways to solve this question, I asked for another approach here.



Whenever a major release of PostgreSQL appears, homebrew will happily install it right away, unless one pins a certain version or installs a specific version from Homebrew/core.


So far so good, but the problem is that PostGIS also keeps updating and eventually will be incompatible with pinned or Homebrew/core versions of PostgreSQL. Again, one can pin PostGIS, yet Homebrew/core does not offer previous versions. From time to time I do brew cleanup and somehow I am now left with a broken PostGIS database. Ok, so "time to jump to PostgreSQL 10.1", I though, but pg_upgrade now gets stuck with the error



ERROR: could not access file "$libdir/postgis-2.3": No such file or directory



That means I can't back nor forth: I cannot figure out how to get PostGIS working with postgresql@9.6 and I cannot pg_upgrade to the new version of PostgreSQL (10.1).


Therefore my question can be answered in two ways. Here I am asking:


How to install PostGIS with older PostgreSQL versions, i.e. postgresql@9.6, using homebrew?


The second approach is described at Upgrading PostgreSQL with broken PostGIS installation ("$libdir/postgis-2.3" error)?





Tried so far:



  • This gist offers a solution and it might actually work. However, I had to use ln -s instead of cp -a. It seems to be working now, but unfortunately I have the problem that my database is already on 2.4 and there seems to be no downgrade path. Will keep trying.




Exporting table from ArcGIS Desktop to Excel file?


Is there a way to export a table with selected records to an Excel spreadsheet?


I was trying to export just to .txt file but ArcGIS gives me an error message.



An error occurred exporting the table




The other problem I have is with coordinates. I want them only with one locus after the comma, I mean like the format: xx xx xx,x not like I have now xx xx xx,xxx.


alt text


The version of my software is 9.3.



Answer



Recent versions of ArcGIS, since 10.2 and we are now at 10.5, have a set of tools named Excel Toolset which has a Table To Excel tool which:



Converts a table to a Microsoft Excel file.



You can convert to excel directly from this tool. It does have some limitations. For example Excel 2010 onward supports 16,384 columns but ArcGIS will only export 256 when I last checked it.


I understand the need to go directly to Excel from a Geodatabase as using the DBF route will truncate your field names. If you are working in a shapefile I am not sure you gain much additional functionality.



See Jay's answer for the Parse command section of the answer as this still stands.


pyqgis - QGIS plugin: Problems importing resources (resources_rc) file - plugin doesn't load - PATH problems?


I'm building qgis plugin and I can't find solution for this error.


File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 478, in _import
mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named resources_napoved_rc

For everybody that will be asking I've built python resources file:


pyrcc4 -o resources_napoved_rc.py resources_napoved.qrc


I still can't find a way to make it work. I always get same error.


On top of the script I have:


import resources_napoved_rc.py


I'm using ui file directly from qtbuilder. Any ideas how to go forward ? I'm assuming this must be some kind of path problem or something similar.




Converting ESRI Raster Dataset to QGIS compatible Dataset?


I am trying to convert an ESRI .gdb "File Geodatabase Raster Dataset" (Format: FGDBR) to an QGIS compatible File with or without using ArcMap.



Answer



You will need to export the image from file geodatabase into a raster file, as I don't believe QGIS can open the file geodatabase raster directly.




  1. In ArcMap Catalog window (or in ArcCatalog itself) browse to your raster dataset in your file geodatabase

  2. Right-click the raster dataset, and select Export > Raster to different format

  3. Select a path to save the raster file to, and choose a raster format that can be opened in QGIS.

  4. Open QGIS and add the new raster file




See GDAL Raster Formats for a list of raster types that should be able to be opened in QGIS. TIF, PNG, JPG are possibly the more common types.


Not all the raster formats in that link are exportable from ArcGIS. These are the formats exportable using the Raster to Different Format tool in ArcMap:


enter image description here



Reading in local GeoJSON file for OpenLayers


I am kinda new to OpenLayers and not (yet) very experienced in JavaScript.


I have a WordPress blog and I want to replicate and adjust the following script to my needs: http://jsfiddle.net/eqsrnjos/1/


With one big exception: I need to load the GeoJSON file from the server instead of defining it in the top (so replace the var geoJson).


However the most straight forward way does not seem to work for me:


function getData() {
return jQuery.getJSON("/json/studyCentroids.json");

}

var features = new ol.format.GeoJSON().readFeatures( getData() );

In Firebug I can see that the GeoJSON file is properly loaded, but it does not seem to appear for some reason?



Answer



var vectorLayer = new ol.layer.Vector({
source: new ol.source.Vector({
format: new ol.format.GeoJSON(),
url: 'http://yourserver.com/yourgeojsonfile.js'

}),
style: new ol.style.Style({
image: new ol.style.Circle( /** @type {olx.style.IconOptions} */ ({
radius: 20,
fill: new ol.style.Fill({
color: '#ffff00'
})
}))
})
});


http://jsfiddle.net/expedio/9hnjok8n/


Thursday 29 June 2017

Understanding difference between Coordinate System and Projection?


Can somebody please explain what is the difference between the Coordinate system (WGS 84 for example) and a Projection (Universal Transverse Mercator for example)?


What is the difference between a projected coordinate system and projected CRS



Answer



Both examples are coordinate systems. The difference is that WGS 84 is a geographic coordinate system, and UTM is a projected coordinate system. Geographic coordinate systems are based on a spheroid and utilize angular units (degrees). Projected coordinate systems are based on a plane (the spheroid projected onto a 2D surface) and utilize linear units (feet, meters, etc.).


More here: Difference between Geographic and Projected coordinate systems?



To answer your second question, a coordinate system (whether geographic or projected) and a coordinate reference system refer to the same thing.


raster - How to de-noise a DSM?


I have created a DEM, DSM, and nDSM from a LiDAR Point Cloud using SAGA-GIS. I extracted the appropriate returns from the point cloud, then converted to grid, and then closed gaps to fill in any nodata regions.


I'm finding that my resultant rasters are still pretty noisy. For example, a region of pixels in the DSM that is obviously a tree or cluster of trees has several pixels with a very similar value to ground instead of the surrounding higher value pixels.


I've tried doing a majority filter on the DSM, but it didn't help.



Any suggestions for other filters that might help create more contiguous regions of pixels?




geodesy - Moving dead ahead: rhumb line, great circle or none


I think I understand the whole rhumb line and great circle idea but I can't seem to find an answer on this question. Suppose I am driving a car (i don't use ship because I don't want to put the constant compass bearing idea on the table) and that I can move for quite a long distance on land (across Russia for example) without any obstacles. If I keep the steer fixed so that the car moves dead ahead will I be following a rhumb line, a great circle or none of these.


I've found this on wikipedia:



If you were to drive a car along a great circle you would hold the steering wheel fixed, but to follow a rhumb line you would have to turn the wheel, turning it more sharply as the poles are approached




but everything else I've found seems to contradict this. And yet I'm not sure if keeping the steer fixed means fixed dead ahead or fixed turned at 30 degrees left for example.


EDIT After a request to give some examples of contradicting information found on the internet I edit my question to quote an example:


http://www.maritimeprofessional.com/Blogs/Maritime-Musings/February-2011/Rhumb-line.aspx



The problem for navigation is that sailing a true great circle would require constant changes in the course steered



I now understand that what the writer of this text probably refers to when writing course steered is holding the compass indication unchanged. But this way of phrasing it is (at least to a not native english speaker like me) very confusing. In my mind "sailing a true great circle would require constant changes in the course steered" means that i have to constantly turn the ship's steering wheel in order to reach my destination.


I had found some similar examples which I can't find again now. But the greatest problem was I couldn't find any other source stating clearly that moving dead ahead equals following a great circle.


However, stating that "everything else I've found seems to contradict this" was a bit of an exaggeration and I apologize for it.




Answer



If the wheels are pointing straight ahead, you are taking the shortest route toward some point directly ahead and thus following a great circle (or geodesic).


If you wish to cross each meridian at exactly the same azimuth each time (and thus follow a rhumb line), you would have to gradually steer slightly more towards the meridian as you progress. Generally, rhumb lines that head between 0° and 90° spiral toward a pole. Parallels of latitude are special cases, with azimuth staying at 90° or 270°.


More special cases to the above are if you happen to be driving straight along a meridian or the equator, which are both great circles and rhumb lines.


Cannot find layer error with Execute SQL in QGIS Graphical Modeler


I'm trying to sort a shapefile by one of its attributes using "Execute SQL", in the QGIS Graphical Modeler.


Using a vector as input is working fine, but when I use a layer generated by an algorithm I have this error message : "virtual: Cannot find layer OUTPUT_shp20180419141642572".


Can someone help me with this ?


To sum-up : here is my problem





arcpy - Using Vector Ruggedness Measure (VRM) tool from ArcScripts at latest version of ArcGIS for Desktop?


I need to use the Vector ruggedness measure (VRM) tool available from http://arcscripts.esri.com/details.asp?dbid=15423. However, I am using ArcGIS 10.2 for Desktop and these tools are made for 9.0, 9.1 and 9.3 only. I have tried to run this tool regardless and get the error:



Traceback (most recent call last): File "C:\Users\jc221340\Documents\project\Processed data\My_tools\Downloaded\ruggedness.py", line 20, in gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx") RuntimeError: Object: Toolbox C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx does not exist


Failed to execute (Ruggedness(VRM)).




How can I use this tool?




javascript - Why the icon can't display with my JSON data on Leaflet



I tried many time, but the icon didn't show on the map. I want to know what the problem on my code.


This is my code:






integrity="sha512-puBpdR0798OZvTTbP4A8Ix/l+A4dHDD0DGqYW6RQ+9jxkRFclaxxQb/SJAWZfWAkuyeQUytO7+7N4QKrDh+drA=="
crossorigin="" />









































Answer



The drawn feature are added to the layer named vector. See draw interactions constructors.


// Get the array of features
var features = vector.getSource().getFeatures();


// Go through this array and get coordinates of their geometry.
features.forEach(function(feature) {
console.log(feature.getGeometry().getCoordinates());
});

QGIS: ESRI shapefile .shp no longer appears in list of file types under Add Vector Layer


I recently installed QGIS 2.18.4 after using ArcMap for many years. I want to add several existing ESRI shapefiles (.shp) created in Arc, to a QGIS project. To add a new layer in QGIS, I select the "Add Vector Layer" button and then click "Browse" to go find the layer. Then I navigate to the proper folder. However, the drop down menu at the lower right of the "Open an OGR Supported Vector Layer" window does not include ESRI Shapefile on the list. The first one on the list is GeoJSON. Everything after "G" appears to be okay (all the way down to "X-Plane/Flightgear") but any file formats from A through F have been cut from the list. Since "ESRI" starts with "E" it does not show up on the list. Strange thing is when I first installed QGIS, ESRI Shapefile appeared on this list. So I must have accidently done something to change how that list is populated and I want to know how to change it back.


If I type in the name of the ESRI shapefile that I know is in the folder, it and all other files starting with those letters appears in the prompt window that appears. And QGIS will correctly load the Shapefile. However, I want to be able to see all those ESRI Shapefiles without having to know and type in the first letter so it will appear in the prompt. I want it to appear in the file type drop down (or drop up?) menu.




georeferencing - Converting longitude and latitude coordinates to Indian Grid system?


How to convert lat/long coordinates to Indian Grid system?


I've googled for a couple of days but didnt find the exact algorithm for this.


and i got a link Indian Grid System.


but i didnt get a picture how the zones(there are 9 zones) are divided among india and how to convert the coordinates?





geoserver - Is there a way to check and repair a GeoPackage GPkg file?


I have downloaded imagery in GeoPackage format from My DigitalGlobe EV which is used offline on a non-internet connected system. However the GPkg seems to be non-standard.


I'm trying to load the GPkg with GeoServer with the GeoPackage extension. However, when I add a GPkg data store I get a Java exception, which I've filed a bug on.


I've tried using GDAL_translate to convert it to a GeoTIFF, and even convert it from GPkg to a new GPkg, but both cause gdal_translate to hang with 100% CPU usage. GDAL also complains about an invalid 'application_id' in the GPkg file.


Is there a script to repair or rebuild a raster GeoPackage file? Before I contact DigitalGlobe, is there anything else I can try?





pyqgis - Selecting and updating attribute values using Python console in QGIS?


I am working with Landusefc and Class columns (ref table).


What I want is: Selecting all water from "Class" and giving a value say '8' in the "Landusefc" column for these selected features.


Table looks like


Test table I have tried the following code, it works till the selection part but doesn't seem to change the attribute value using the ChangeAttributeValue.


from qgis.core import *
import processing

layer=processing.getObject('test1')
query= ' "class" = \'Water\' '
selection = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))
layer.setSelectedFeatures([k.id() for k in selection])

#using change Attribute to change the value of the selected feature

layer.startEditing()
for feat in selection:
layer.changeAttributeValue(feat.id(), 5, 8)


layer.commitChanges()

Answer



SOLUTION 1:


You have to remove layer.setSelectedFeatures([k.id() for k in selection]).


Because selection = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query)) line gives you a QgsFeatureIterator, that's an iterator. When you use that iterator once, you reach the end of the iterator. Then, if you use it again, you get nothing.


In your case, when you use selection variable (it's an iterator now) in [k.id() for k in selection], you reach the end. If you use selection in for loop, since you are at the end of the iterator, you don't have any feature now. If you want to be certain, you can add a print("Foo") line to for loop.


SOLUTION 2:


Change for feat in selection: into for feat in layer.selectedFeatures():. This solution works only on QGIS 2.x.


Sunday 25 June 2017

How do I verify a file is a GeoPDF?


I received some files from a third party that are supposed to be geoPDF. Is there a way to open a file (programmatically with .NET or Python, with a text editor, etc.) and verify it was created as a geoPDF? Does a geoPDF have bytes or an internal structure we can check?




qgis - How to add my own GRASS algorithms to Processing?



I need to use algorithm r.bioclimatic from GRASS in QGIS but I see that this algorithm don't exist in Processing Toolbox. A search a way of add this algorithm to Processing and in GRASS and Sextante page if find this way:


HOWTO adding own GRASS algorithms to Sextante (simply by adding a new description file)


I create a description file according to the rules but now I don't know how to run and if this will run in GRASS.


The first parameters of r.bioclimatic (tmin, tmax....) are columns of tables? I don't know which type of input I have.



Answer



r.bioclim (I suppose, because r.bioclimatic don't exist), is an GRASS GIS 7 addon


If you want to create or modify an algorithm, read the file ../processing/algs/grass7/grass7.txt that begins with "A short guide for creating and editing GRASS GIS 7 algorithms"


But, how do these commands work ?



  • to use the GRASS7 commands in processing,you need to first specify the GRASS GIS installation (in Processing/options) and processing will look for scripts that are in the folders bin and scripts of your GRASS installation (GISBASE environment variable) in the script ../processing/algs/grass7/Grass7AlgorithmProvider.py


  • the script is executed with the Python module subprocess in the ../processing/algs/grass7/Grass7Utils.py script


The big problem is that the addons are not installed in these folders (GRASS_ADDON_PATH environment variable) therefore you can simply use them with processing


Conclusion ? QGIS is not GRASS GIS and if you want to use a "specialized" algorithm (addon), use GRASS GIS


qgis - Projection system differs


I want to define projection for a Shape file. The projection system is: "NAD_1983_StatePlane_Montana_FIPS_2500" which is referred in Qgis 1.9 as EPSG:102700.



The same projection system is not available in Qgis 1.8. Here it is designated as "NAD 83/Mountana(ft)" and has EPSG:2256 and "NAD83/Mountana" and has EPSG:32100.


My problem is if I want to define projection in Qgis 1.8 what projection system should I use? EPSG:2256 or EPSG:32100?



Answer



2256 and 32100 are the same except 32100 is in meters. I looked up NAD_1983_StatePlane_Montana_FIPS_2500 and the proj4 is the same as 2256 except that lat1 and lat2 are switched, which I don't believe really matters.


So 32100 should work. Otherwise, create a prj file for the shapefile with the following WKT:


PROJCS["NAD_1983_HARN_StatePlane_Montana_FIPS_2500",GEOGCS["GCS_North_American_1983_HARN",DATUM["NAD83_High_Accuracy_Regional_Network",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["False_Easting",600000],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-109.5],PARAMETER["Standard_Parallel_1",45],PARAMETER["Standard_Parallel_2",49],PARAMETER["Latitude_Of_Origin",44.25],UNIT["Meter",1],AUTHORITY["EPSG","102300"]]


I got this from http://www.spatialreference.org/ which is a great reference for projection problems.


data - Where can I find Road shapefile/e00/etc for the entire world?


I am hoping that there is a single source for this data, ie. 1 geodatabase. If not, can someone point me towards the sources that would make up the entire world's road transportation network?


I have come across references to the ESRI Digital Chart of the World, but I can't find any place to download it.


Ideas?



Answer



Digital Chart of the World is quite old now:


The correct location is 'http://ortelius.maproom.psu.edu/dcw/' (but unavailable currently)


Other locations are available but not in e00 or shapefile http://statisk.umb.no/ikf/gis/dcw/#DCW


But you might be better getting the Open Street Map data for Roads. But will be large amounts of data for the world.



"This is almost 20 GB compressed or 150 GB uncompressed XML."


http://wiki.openstreetmap.org/wiki/Downloading_data#All_data_at_once


For Shapefiles you best bet is http://download.geofabrik.de/osm/


arcgis desktop - Most Effective Format to Manage Aerial Photography for Viewing Only


What is the most effective format to manage aerial photography for viewing in ArcGIS?


I manage several datasets of Aerial Photography at a few different companies. I have backups of the original photography, so I'm not worried about keeping the original format in any way. What I need to do is create a small dataset, as too not take up much space on the server, that also draws as quickly as possible in ArcGIS.


So far I had been using GRID's in a File GeoDatabase, but that was just because I has assumed it was the standard. I then switched to TIFF's which were smaller, but seemed to draw slower. Before I do another mass conversion, I'm hoping someone can offer guidance on this. I've heard some great things about JP2's and that's where I am likely heading next.


Also I don't want to reduce cell size, but other than that I have no restrictions in any answer.



Answer



Generally I'd go with ECW over JP2 as the refresh rates for the ECW are better. There is MrSID to consider, but it is expensive - possibly more than ECW provided you're appropriately licensed, but at the extreme compression level (95% to 99%) they are clearer than similar ECW.


JP2 offers open source compressors, lossless compression and 4 band support (RGBA). A bit of a loss with aerial photography shouldn't matter too much; The alpha band can be great if you've got two scenes with irregular edges that need to be superimposed.



If you build the pyramids (use GDALAddO to build ovr files) for TIFF files they can be very fast too, and when you zoom in it depends on your compression (LZW, Zip, Deflate... many supported) as to how fast they are. There is a trade off between smaller sizes and CPU usage by using more intense compression and there is a possibility that some software will not like a particular compression method.


In the end JP2 is the best free format for aerial photography and IMHO ECW is the best paid format.


Saturday 24 June 2017

python - Getting the projection of a LAS file using liblas?


I am using the python api of liblas to read/write las files. I would like to get the spatial reference information of the las file using liblas, but have been unable to figure it out.


Can anyone offer some advice on how to go about getting this information?



Answer



To be honest, I think the easiest way would be to shell out to lasinfo with the --xml argument, and then use ElementTree or some such to pluck out what you need from the XML.


You can reach all this stuff from the Python bindings too, but it's a bit of a mess. In short, open the file, fetch the header, and fetch the srs, and then ask for the wkt of the srs. This all assumes you have libLAS + GeoTIFF + GDAL enabled, and degrades poorly if you don't have all three of those things linked together.


Additionally, if all you want to do is interact with LAS files, laspy is a much better alternative to libLAS at this time. laspy doesn't allow you do anything with the coordinate system stuff in LAS other than to fetch the GeoTIFF bytes out of the file, but it is a full-featured library in every other way. It's pure Python + numpy, easy to deploy, and behaves naturally, unlike libLAS (disclaimer, I'm the libLAS author and I helped bootstrap laspy as an effort to obsolete libLAS).


Hope this helps,


Howard


arcgis 10.0 - Why are Character Marker Symbol Dialog not showing up in the Symbol Property Editor dialog?


I just installed ArcGIS 10.2.1 on a brand new machine and for some reason Character Marker Symbols are not loading into the Symbol Property Editor.


Is anyone experiencing the same issue? Is there a solution for this?


enter image description here



Answer



I apologize for ansering my own question this quick but I checked on the ESRI forums and a similar problem is documented there: http://forums.arcgis.com/threads/9869-symbol-property-editor-will-not-display-character-marker-symbols


The culprit seems to be ATI video card with their application called Hydravision. I have an AMD FirePro W7000 (FireGL V) but it also does have HydraVision.


The solution is to disable the disable the "enable dialog repositioning" option.


arcgis desktop - The updateDisplay action in Agent Analyst does not update my ArcMap display


I use ArcGIS 10.4. and Win 8.0. .I have a problem with Agent Analyst in a model-level action. I did all steps of the Exercise 2f in the book.The model start running, the Repast output suddenly stopped working, and an error shows " refresh has stopped working", "check online for the solution...".The updateDisplay action in Agent Analyst does not update my ArcMap display. When I deleted the agent's updateDisplay, I didn't receive any Error.The attached file is the print screen of the Error.


enter image description here



Answer



I created a virtual C: Drive with Windows 7 on my Windows 8 and installed ArcGIS 10.2 on virtual Drive. then I went through the structures mentioned in this Agent Analyst tool interface will not open on ArcGIS for Desktop? then the problem fixed.


Friday 23 June 2017

qgis - Finding polygons without right angles using Open Source GIS or ArcGIS for Desktop?


Now we are digitizing some buildings in specified area.


The obligatory rule for this work - in most cases buildings should have right angles.


We are using QGIS with CAD tools for this work but sometimes we make mistakes and create polygons with an irregular shape.


Does anybody know how can we find such polygons without right angles using open source GIS or ArcGIS?



Answer



I don't know of an existing tool to do this but you can write one in ArcPy or using GDAL/OGR along the following lines:



  • For each polygon...


    • Get the geometry

    • Follow the winding of the polygon and calculate the internal angle at each vertex

    • If any angle is not 90 degrees, add the FID/OID (or some other attribute) to a list of rejects



  • Print the list of rejects


qgis - Display coordinates of points in DMS format in print composer


In QGIS 2.14.2 I have a series of polylines of which I want the start and endpoint to be labeled with its coordinates in WGS84 (EPSG 4326) like this: 51*15'23.40"N, 5*17'56.58"E


I thought I figured a workaround by:



  1. Creating new point shapefile and creating points manually at start and endpoint of each polyline.


  2. Open field calculator and create a virtual field "X" containing $x and another containing $Y (Decimals)

  3. Label them by X||'\n'||Y


However the field calculator gives me an output like this X:5.26344964931 Y:51.20551633106


How do I display them in DMS format?



Answer



You're in the right direction. What you need to do is use the QGIS expression builder to concatenate a string with the degrees, minutes, seconds, and any symbols you want in between.


This example calculates the $x value, and will work for positive or negative values.


(CASE WHEN $x < 0 THEN '-' ELSE '' END) || floor (abs($x)) || '° ' || floor(((abs($x)) - floor (abs($x))) * 60) ||'\'' || substr( (tostring((((abs($x)) - floor (abs($x))) * 60) - floor(((abs($x)) - floor (abs($x))) * 60)) * 60),1,5) || '"'


You can change the symbol signs, like °, ' and '' to fit any style you need.


Do note though that the resulting field has to be a string.


Seeking Esri tools to design geodatabase?



There is a confusing multitude of tools written by Esri and others (some supported, some unsupported) for designing and diagramming a geodatabase.


Here is an attempt at creating a definitive list of where to find them.


Related to What are convenient uml tools to create geodatabase models?



Answer








  • Esri CASE tools with Microsoft Visio, CASE Support. Supported by Esri.

  • ArcGIS Geodatabase Design with UML (for 10.0). Add-in to Enterprise Architect (third party). Unsupported by Esri.

  • ArcGIS Diagrammer (for 9.x, 10.0, 10.1, 10.2). Visual editor (not uml). Schema and data reporting. Unsupported by Esri.

  • Geodatabase Diagrammer (for 9.x, 10.x) Requires Visio. Originally created by Mike Zelier. Unsupported by Esri.

  • X-Ray for ArcCatalog (for 10.0, 10.1, 10.2) Created by Steve Grise of Vertex3. Unsupported by Esri.

  • Geometric Network Configuration Manager (for 10.2) Save, edit and restore geometric network configuration. Unsupported by Esri.


Thursday 22 June 2017

arcgis desktop - Loop in ModelBuilder for compositing bands of multiple Images


I'm attempting to automate the "Composite Bands" tool somehow in ModelBuilder, but haven't been able to figure out an effective method.
I have landsat 7 data that is separated into folders and each holds 7 bands (tif). Basically I want to generate a model or script that composites all bands in the raster output for each location, and then saves the result in a output folder. I've prepared a ModelBuilder that take each band and save it in its correspondent folder. However, the input raster just read 1 band, not all 6 bands. Any Solutions, maybe Python? Model Builder


Composite Bands which just i raster uploaded :(



Answer



You mentioned that you are open to a Python solution, which is good because automating this process with Python is much easier and more flexible than a ModelBuilder approach.


First, import the necessary modules


import arcpy, os

Define the workspace that contains all of the folders with the Landsat imagery



arcpy.env.workspace = r'C:\imagery'

Specify where you would like the output to go


outws = r'C:\temp'

List all of the workspaces in the previously defined workspace


folders = arcpy.ListWorkspaces()

Iterate through this list of workspaces and create a new list within each iteration of unstacked raster bands e.g. ['LC80260272014159LGN00_B1', 'LC80260272014159LGN00_B2',...]


for folder in folders:

arcpy.env.workspace = folder
rasters = arcpy.ListRasters("*.tif")
name = os.path.join(outws, rasters[1].split("_")[0] + ".tif")
arcpy.CompositeBands_management(rasters, name)

Landsat files follows this form: LC80260272014159LGN00_B1.tif, so we need to strip off anything after the "_" and use the first basename as the output name. You can do this with various slicing methods and string manipulation in Python.


name = os.path.join(outws, rasters[1].split("_")[0] + ".tif")



The complete script:



import arcpy, os

arcpy.env.workspace = r'C:\imagery'
outws = r'C:\temp'

# list all folders in a directory

folders = arcpy.ListWorkspaces()

for folder in folders:

arcpy.env.workspace = folder
rasters = arcpy.ListRasters("*.tif")
name = os.path.join(outws, rasters[1].split("_")[0] + ".tif")
arcpy.CompositeBands_management(rasters, name)

print "Processing complete"

arcgis desktop - How to show coordinate on the polygon



how to it


I try to set and go to search by help GIS but I couldn't find an answer


How to show coordinate on the polygon but I don't show it on dataframe




qgis - Advanced colour palette for raster data


I'm working with satellite date namely chlorophyll and sst. Each geophysical parameter generally has standard colours to represent them. I'm importing geotiffs into Qgis and I'm having difficulty finding the colours needed and a way of styling the raster data.


Is there a way to add colour maps when in the raster properties?


Can I show a colour bar in print composer or in the legent for these colours.


I found a way of changing each colour individually but I'm hoping to find a standard range of colours like the pseudocolour palette.


Any help is greatly appreciated.


Regards




Wednesday 21 June 2017

accuracy - How accurate are measurements in Google Earth?


Can anyone point to a resource that lists the accuracy of measurements taken in Google Earth?



Answer



Apparently, the measurment ruler is not accurate over long distances > 5,000 km.


Google's offical stance says, "makes no claims as to the accuracy of the coordinates in Google Earth. These are provided for entertainment only and should not be used for any navigational or other purpose requiring any accuracy whatsoever".


QGIS - adding timestamps to spatialite data during data collection


I am trying to do field data collection using QGIS. One of the features that I regard as necessary for data management is timestamps. For example, a timestamp can indicate that an inspection record has been updated.


I have figured out how to add triggers using Spatialite to add a timestamp on update, however the timestamp is added when the edited records are saved and not when the identify dialog box is closed. This means that timestamps can be out by many minutes. One solution would be to force the updated record to be saved by using a custom form with an OK button (Can do the form with Qt-Designer but don't know how to add the save command to the button yet.) However my preferred solution would be to take a timestamp for new records, when the point or line is created. The reason for this is if I am surveying a road, I can stand on the road, collect the point and them move to the side to complete data collection.


I have been unable to work out how to get Spatialite to run a trigger that puts in the start of feature editing. Perhaps this is also something that needs to be done by code in a custom form. Is there an On-dirty event or signal in QGIS that indicates when someone starts editing a record and that can be used to capture the time? Could I also ask what the Python code (and what to hook it to) that would put such a timestamp in the timestamp field.


AndrewM



Answer



Note: This will only work in 1.9 (2.0) due to a bug in 1.8 not passing the full feature correctly


Take this code:



from PyQt4.QtCore import QDateTime

def onopen(dialog, layer, feature):
if feature.id() < 0:
# This is a new record because it
# doesn't have a vaild id yet
feature['createdon'] = QDateTime.currentDateTime()

Save it somewhere QGIS Python can find it, normally in the same folder as the project is good. Lets call it datestamper.py


In Layer Properties -> Fields set the Python Init Function to datestamper.onopen



enter image description here


Now when ever the form is opened, auto generated or custom, that function will be called and the 'createdon' field will be set to the current date time.


data - Seeking Sentinel 2 Images for Nigeria




I need Sentinel 2 imagery of Nigeria, I signed up in Copernicus Scihub but typing Nigeria into the text box no results appear. I must be doing something wrong, how can be possible that they do not have images for Nigeria?




installation - QGIS 2.14 fails to install on Ubuntu 14.04


My previously working installation of QGIS 2.14 on Ubuntu 14.04 suddenly failed to start QGIS. I have been trying to reinstall QGIS without success. I have tried several ways to completely uninstall first, including this -


sudo apt-get autoremove
sudo apt-get autoclean
sudo apt-get -f install
sudo apt-get autoremove qgis
sudo apt-get --purge remove qgis


Then, on trying to reinstall QGIS, install fails with this error -


Setting up qgis-providers (1:2.14.0+20trusty-ubuntugis) ...
/usr/lib/qgis/crssync: symbol lookup error: /usr/lib/libqgis_core.so.2.14.0: undefined symbol: initGEOS_rdpkg: error processing package qgis-providers (--configure): subprocess installed post-installation script returned error exit status 127

So qgis-providers is not getting configured. I have no idea what the undefined symbol issue might be.


I was using the qgis debian repository (http://qgis.org/debian) + the Ubuntugis repository (ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu). Also tried qgis ubuntu repos (http://qgis.org/ubuntugis). Tried only debian qgis repos (with ubuntugis turned off. All combinations fail configuring qgis-providers with same error.


I was able to install qgis 2.0 from the standard Ubuntu repository apparently without problem, but I find that my QGIS project files saved from 2.14 will not open in 2.0.


This has me stumped. I thought for sure I could just reinstall to recover, but after several frustrating hours, I could sure use some help or clues to some resolution.




Custom functions for SAVE and CANCEL on Leaflet.draw edit control


I have added an edit control on my map which is working fine (On the Frontend);


var editControls = new L.Control.Draw({
draw: false,
edit: {
featureGroup: this.updateGeomLyr
}
}).addTo(this.map);


enter image description here


Now, when I save or cancel my changes, the event "draw:editstop" is called


this.map.on("draw:editstop", e => {
this.updateGeomLyr.clearLayers();
this.overlapLyr.clearLayers();

//if saved;

//if cancelled,

});

But I want separate backend functions to be called on SAVE ad CANCEL. I couldn't find any attribute in the event "e" which shows if save is called or cancel.




Calculate scale of a raster


A raster has a cell size of 200 meters. What is the corresponding scale of the raster?



New to GIS, using ESRI, semester is almost over but this one question is confusing. I didn't think you could calculate the scale of the raster in that you are given cell size of 200 meters and to my knowledge, or lack there of, the corresponding scale would be whatever the data set is using or that you adjust the scale to.




python - Dynamic Text in ArcMap 10.0 and arcpy - getting user name in specific format


I need to put in a specific text into a text element. The initials of the user who built the map.


With arcpy I am trying to set the text element with dynamic text. If I put



 

in the user text element it gives me, correctly, the user's user name. I want to use the user name in a python dictionary to put the users initials instead.


dictUserElmt = { "jperez" : "JP" }


Using arcpy to read the text of my text element, textElement.text, returns (again correctly)


u''.  

How do I get the actual output of the dynamic text "jperez" so I can use my dictionary to have the text element set to "JP".


I'm sure my wording is confusing because I am confused on how to approach. Thanks.



Answer




Is something like below what you are looking to do? For a given username look up and return their initials from the dictionary?


>>> d = {"jperez":"JP", "cooperc":"CC"}
>>> user = "jperez"
>>> if user in d.keys():
... inits = d[user]
...
>>> inits
'JP'

Now, that said, are you getting the author name from the MXD properties?



>>> import arcpy
>>> mxd = arcpy.mapping.MapDocument("CURRENT")
>>> a = mxd.author
>>> print a
>>> a
u''
>>> a = mxd.author
>>> print a
Chad Cooper


If so, that has to be set by the user (you can see I got u'' on my trial before I explicitly set my name as the MXD author) and saved in the MXD I believe. In your code, you have your name as jperez, which looks like a Windows login maybe? I wonder how using the getpass module to get the current login name would work? That way you could go by those, which pretty much stay the same.


>>> import getpass
>>> getpass.getuser()
'chad'
>>>

Tuesday 20 June 2017

Adding an external WMS service using REST service with Geoserver


I'm trying to add an External WMS layer to geoserver with a PHP page, using the REST service of Geoserver. I've already read the tutorials for creating the geoserver features (http://docs.geoserver.org/stable/en/user/restconfig/rest-config-api.html), but i couldn't find the way to create a new store and layer with an external WMS service.


in other words: I have the URL of the external WMS service (=store), the layer i want to publish (=layer) and the epgs i want to use.


Is there a way to load all of these information in the geoserver schema?




Finding middle point (midpoint) of line in QGIS?


How can I locate the coordinates of the middle point (midpoint) of a line feature in QGIS?




Adding TauDEM provider to QGIS 3?


Is there a way to add TauDEM to the latest QGIS 3.0?


I can't see it in the providers tab as it was in my previous 2.18.14 version so I am unable to activate it. I use it pretty regularly so I would really like to get it added to this latest release as soon as possible to test out its new features.



Answer



TauDEM was removed from QGIS 3.0 as you mentioned in the changelog of QGIS 3.0: remove TauDEM provider from core Processing. The only solution for you for the time being is to use QGIS 2.18.17 (New LTR) which still supports TauDEM.


I am not sure if TauDEM will be added in the future or not, but at least you can still use it with QGIS 2.18 version.


web mapping - Maximum number of point features in an OpenLayers vector layer


In your experience, how many point features can be added to an OpenLayers vector layer (new OpenLayers.Layer.Vector("Point Layer")) before it goes unusably slow?


My use case is to display points from a database table. The user can decide which time frame to visualize. Therefore the result can be from very few to potentially 100,000s of points. I'd like to introduce a reasonable limit and warn the user if his query would return more features.



Answer



I don't have a definitive answer for you but you I put together a page where you can play around with different numbers of points on an OL map: http://derekswingley.com/lab/olpts/


postgis - Extract lat, lon Points from Polygon



I have PostGIS database, and want to extract lat, lon pairs sequence from Polygons (and MultiPolygons). As a very fresh PostGIS user, I'm able to just extract the Polygons:


SELECT st_astext("GEOMETRY"::geometry) FROM my_postgisdb;


which outputs Polygons and MultiPoligons.


How should I extract x, y pairs from these structures?



Answer



You can use ST_DumpPoints to get the points and then ST_AsLatLonText to get nicely readable coordinate output. It doesn't care if the geometries are multi* or not.


Alternatively, you could use ST_Boundary to get the polygon vertices, but you'd need to use ST_Dump first for this to work for multipolygons too.


spatial statistics - Determining if trees within forest gaps are clustered using R?


The attached dataset shows approximately 6000 saplings in approximately 50 variable sized forest gaps. I am interested in learning how these saplings are growing within their respective gaps (i.e. clustered, random, dispersed). As you know, a traditional approach would be to run Global Moran's I. However, aggregations of trees within aggregations of gaps seems to be an inappropriate use of Moran's I. I ran some test statistics with Moran's I using a threshold distance of 50 meters, which produced nonsensical results (i.e. p-value = 0.0000000...). The interaction among the gap aggregations are likely producing these results. I have considered creating a script to loop through individual canopy gaps and determine the clustering within each gap, although displaying these results to the public would be problematic.


What is the best approach to quantify clustering within clusters?



enter image description here



Answer



You do not have a uniform random field, so attempting to analyze all of your data at once will violate the assumptions of any statistic you choose to throw at the problem. It is unclear from your post if your data is a marked point process (i.e, diameter or height associated with each tree location). If this data is not representing a marked point process I have no idea how you applied a Moran's-I. If the data only represents spatial locations, I would recommend using a Ripley's-K with the Besag-L transformation to standardize the null expectation at zero. This allows for a multiscale assessment of clustering. If your data has an associated value, then your best option is a local Moran's-I (LISA). I would actually look at it with both statistics. Regardless of your choice, you will still need to loop through each individual site to produce valid results. Here is some example R code for a Monte Carlo simulation of Ripley's-K/Besag's-L using the built-in redwood sapling dataset. It should be fairly straightforward to modify this to loop through your sites and produce a graph for each one.


# ADD REQUIRED PACKAGES
require(sp)
require(spatstat)
options(scipen=5)

# USE REDWOOD SAPLING DATASET
spp <- SpatialPoints(coords(redwood))


###################################################
###### START BESAG'S-L MONTE CARLO ANALYSUS ######
###################################################
# CREATE CONVEX HULL FOR ANALYSIS WINDOW
W=ripras(coordinates(spp))

# COERCE TO spatstat ppp OBJECT
spp.ppp=as.ppp(coordinates(spp), W)
plot(spp.ppp)


# ESTIMATE BANDWIDTH
area <- area.owin(W)
lambda <- spp.ppp$n/area
ripley <- min(diff(W$xrange), diff(W$yrange))/4
rlarge <- sqrt(1000/(pi * lambda))
rmax <- min(rlarge, ripley)
bw <- seq(0, rmax, by=rmax/10)

# CALCULATE PERMUTED CROSS-K AND PLOT RESULTS

Lenv <- envelope(spp.ppp, fun="Kest", r=bw, i="1", j="2", nsim=99, nrank=5,
transform=expression(sqrt(./pi)-bw), global=TRUE)
plot(Lenv, main="Besag's-L", xlab="Distance", ylab="L(r)", legend=F, col=c("white","black","grey","grey"),
lty=c(1,2,2,2), lwd=c(2,1,1,1) )
polygon( c(Lenv$r, rev(Lenv$r)), c(Lenv$lo, rev(Lenv$hi)), col="lightgrey", border="grey")
lines(supsmu(bw, Lenv$obs), lwd=2)
lines(bw, Lenv$theo, lwd=1, lty=2)
legend("topleft", c(expression(hat(L)(r)), "Simulation Envelope", "theo"), pch=c(-32,22),
col=c("black","grey"), lty=c(1,0,2), lwd=c(2,0,2), pt.bg=c("white","grey"))

installation - Where can I get ArcView 1.0?



The version is correct 1.0. Not 10.0, 1.0. Will it run on Windows 7?



Answer



From the The Surveying, Mapping & GIS blog post blast from the past. (direct download link). You'll need Windows 3.1 or Windows for Workgroups, though I wouldn't be surprised if FreeDOS or Wine would work as well (or better!).


qgis - Toggle Map Tips with PyQGIS?


I want to display "map tips" in QGIS with python for my custom map canvas QgsMapCanvas().


How can I do this? I can switch manually with View/Map Tips but it won't work with my custom canvas in my plugin, it works only with iface.mapCanvas().



Answer



Not sure if there's another way, but this is how I've implemented it.


You need these three functions (make sure you import the required QGIS and Qt4 classes before):


def createMapTips( self ):
""" Create MapTips on the map """
self.timerMapTips = QTimer( self.canvas )

self.mapTip = QgsMapTip()
self.connect( self.canvas, SIGNAL( "xyCoordinates(const QgsPoint&)" ),
self.mapTipXYChanged )
self.connect( self.timerMapTips, SIGNAL( "timeout()" ),
self.showMapTip )

def mapTipXYChanged( self, p ):
""" SLOT. Initialize the Timer to show MapTips on the map """
if self.canvas.underMouse(): # Only if mouse is over the map
# Here you could check if your custom MapTips button is active or sth

self.lastMapPosition = QgsPoint( p.x(), p.y() )
self.mapTip.clear( self.canvas )
self.timerMapTips.start( 750 ) # time in milliseconds

def showMapTip( self ):
""" SLOT. Show MapTips on the map """
self.timerMapTips.stop()

if self.canvas.underMouse():
# Here you could check if your custom MapTips button is active or sth

pointQgs = self.lastMapPosition
pointQt = self.canvas.mouseLastXY()
self.mapTip.showMapTip( self.layer, pointQgs, pointQt,
self.canvas )

You would also need to set a layer, and probably a display field:


self.layer = layer # For MapTips
self.layer.setDisplayField('NMG')

Finally, initialize the createMapTips function, for example, after you set your base layer and its display field.



self.createMapTips()

You should be able to see map tips on your custom mapCanvas like this:


enter image description here


Notice that you could need some GUI controls to help users use and configure map tips. For example, a button to enable/disable them and a comboBox for users to choose the display field, but, of course, it would depend on your use case.


spatial analyst - Error: Intersect 2 polygons in Python. Record's geometry not match: 'GeometryCollection' != 'Polygon'



I'd tried replicating intersecting two shapefiles from Python or command line by anna-magdalena which requires intersecting two shapes. I used @gene answer for that:


import os
from shapely.geometry import mapping, shape
import fiona

os.chdir('/media/dir')
# Create schema of values that new shape will contain
schema = {'geometry': 'Polygon','properties': {'area': 'float:13.3', 'id1': 'str', 'id3': 'str'}}

with fiona.open('13fiona.shp', 'w',driver='ESRI Shapefile', schema=schema) as output:

for c1 in fiona.open('1.shp'):
for c2 in fiona.open('3.shp'):
if shape(c1['geometry']c1).intersects(shape(c2['geometry'])):
area = shape(c1['geometry']).intersection(shape(c2['geometry'])).area
prop = {'area': area, 'id1' : c1['id'],'id3': c2['id']}
geom = mapping(shape(c1['geometry']).intersection(shape(c2['geometry'])))
output.write({'geometry':geom,'properties': prop})

And I got this error:


ValueError: Record's geometry type does not match collection schema's geometry type: 'GeometryCollection' != 'Polygon'


After check the structure of objects c1, c2 and geom I got this:


Structure of c1 polygon:


{'geometry': {'type': 'Polygon', 'coordinates': [[(-74.71213794799996, 1.2098971930000744), (-74.71 ...

and c2 polygon:


{'geometry': {'type': 'Polygon', 'coordinates': [[(-74.69988807399994, 1.2222...

I found that intersection polygon assigned to geom objetc have this structure:


{'type': 'GeometryCollection', 'geometries': [{'type': 'Point', 'coordinates': (-74.71039521399996, 1.2098918560000698)}, {'type': 'Point', 'coordinates': (-74.70929715499994, 1.2102729000000636)}, ...


The c1 and c2 shapes cand be found here I'm using shapely 1.6b2 version


Any advice?




arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...