On the http://spatialreference.org, for almost all SRS are definied WGS84 Bounds and Projected Bounds.
How to obtain this bounds from postgis? Spatial_ref_sys does not contains this information.
Answer
You'll have to download the EPSG database and pull the information you want from it. They do provide PostgreSQL-compatible downloads. http://www.epsg.org/DownloadDataset
Is it possible to postprocess Trimble .ssf files without Trimble Pathfinder Office (differential correction - our sporadic use does not defend realtime correction subscription) ? We only need the correction utility but the pathfinder package contains a lot more and is too expensive just for this. I have not yet found open source software which supports ssf.
There is probably already a post for my problem, but after long searching, I could not find a solution for my issue.
I am working with QGIS und I am using the following projection/datum: MGI / Austria GK East (EPSG:31256) (+proj=tmerc +lat_0=0 +lon_0=16.33333333333333 +k=1 +x_0=0 +y_0=-5000000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs)
I just want to add a layer form the openlayers plugin. For example Google satelite. the Google satelite uses the following projection/Datum: WGS 84 / Pseudo Mercator (EPSG:3857) (+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 +wktext +no_defs).
The Problem that I am facing with, is that I can not get the Google layer in the same projection with my other layers (EPSG:31256)! They are not overlapping correct!
Is there a way to do this?
I'm using pgrouting DLL 1.05 with postgresql 9.1. All required function installed.
Astar and Dijkstra work fine with sample data but "SHOOTING STAR" is not considering "reverse_cost".
Also if i try long route then it gives server disconnected error.
"Do you want to attempt to reconnect to the database?"
And then I need to restart the postgresql service.
I have just upgraded to ArcGIS 10.1 from version 10.0. I am trying to add the SEXTANTE tools to ArcToolbox. I have followed all the instructions from the SEXTANTE website. I have also tried the instructions here and here.
No matter what I do, I keep getting little red crosses when I try to run any of the tools in SEXTANTE.
Note, following the instructions in the above links worked for me in version 10.0 but I just can't seem to get it to work in 10.1.
Does anyone know what might be going on? Are other people having this issue?
Answer
I've written to the SEXTANTE team to ask about this issue. This is the response I got:
Other people have reported that those "tricks" that worked for 10 do not work for 10.1. We haven't tested it, but there are plans for doing some development on the ArcGIS version in collaboration with a team from an american university, and definitely that's our priority, to get rid of that problem when installing.
I'll update this answer when they've fixed the issue.
UPDATE 15 Feb 2013:
I just checked this page again to see if there has been any update and I get a 404 error. Looks like they've pulled support for ArcGIS for the time being.
Ok, likely my bad for not making it to WhereCamp5280 and asking the FortiusOne geeks directly, but what's the process for creating GeoJSON 'tiles' for use in Polymaps?
The Population example over at polymaps.org notes that the data for the demo runs on Google AppEngine, and alludes to "rolling your own" vector tiles, but I can't seem to find any more information...
Anyone got some insight they can share?
Thanks!
Answer
TileStache will definitely do it, specifically the PostGeoJSON provider in the extra providers collection. The main difference between this and what we did for the Polymaps examples is polygon clipping - after seeing how large Alaska can get at higher zoom levels, we clipped all the features in the Polymaps examples to cut down on load times and so on. We also hosted them out of AppSpot to make it possible to add the Access-Control-Allow-Origin header for cross-domain permissions.
PostGeoJSON doesn't clip out of the box, but as it says in the docs I'd be happy to develop the code further if there's interest!
Update: It appears that as of Tilestache 1.9.0 the Vector Provider appears to be favored over the PostGeoJSON Provider.
I have an attribute field in my file geodatabase (GDB) where I'm trying to select all records that contain an underbar in the string field. The issue is that the underbar is considered a SQL wildcard on its own.
To illustrate: I need all data similar to "ID" = '42W_42_42E' or "ID" = '31A_31'
Using "ID" LIKE '%_%' returns all records, which isn't helpful. The perfect solution would simply return a count on the number of records with the underbar in their ID attribute.
Answer
From the page on SQL reference for query expressions used in ArcGIS
To include the percent symbol or underscore in your search string, use the ESCAPE keyword to designate another character as the escape character, which in turn indicates that a real percent sign or underscore immediately follows.
For example, this expression returns any string containing 10%, such as 10% DISCOUNT or A10%:
"AMOUNT" LIKE '%10$%%' ESCAPE '$'
For your case, try the following:
ID LIKE '%$_%' ESCAPE '$'
I'm executing a long process WPS of around 180sec, and when it start after 60sec the WPS (response) is abandoned, and the process continue to run until the end (checked on the log file). I tried the solution proposed by Pywps wiki of reducing the timeout of the Apache server (http://wiki.rsg.pml.ac.uk/pywps/Async_issue) without success !
I found the response which seems to be the only issue : http://lists.wald.intevation.org/pipermail/pywps-devel/2013-April/001598.html But I didn't find how can I insert my "call process" by self.cmd("process.sh...).
Answer
The problem came from the Qgis WPS client, he abandoned the application process lasting more than 60 minutes! I used an openlayers client, it worked well in asynchronous mode, of course, once the timeout instance of PyWPS and Apache have been set up.
I have tried to apply a cloud mask for a landsat 8 collection using code1, code2, and code3 but none will work.
code1 encounters the error of not working even after changing
var mask = img.select(['cfmask']).neq(4) to
var mask = img.select(['pizel_qa_bands']).neq(3)
code2 encounters the problem of not being a time series dataset. I want to be able to extract pixel values for each image over a period of time but this code only gives one output of pixel values.
and I can't adapt my script to code3.
Here is my script and all i want to do is remove the clouds and cloud cover from a Landsat8 image collection
var l8_mayo = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterDate('2000-01-01', '2018-12-31')
.filter(ee.Filter.eq('WRS_PATH', 208))
.filter(ee.Filter.eq('WRS_ROW', 23));
var visParams = {bands: ['B4', 'B3', 'B2'], max: 0.3};
Map.setCenter(-9, 53, 8);
Map.addLayer(l8_mayo, visParams, 'l8_mayo collection');`
Answer
code2 works. Adapt it properly with right parameters (Always check metadata):
var l8_mayo = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterDate('2000-01-01', '2018-12-31')
.filter(ee.Filter.eq('WRS_PATH', 208))
.filter(ee.Filter.eq('WRS_ROW', 23));
var visParams = {bands: ['B4', 'B3', 'B2'], max: 0.3};
Map.setCenter(-9, 53, 8);
Map.addLayer(l8_mayo, visParams, 'l8_mayo collection');
var getQABits = function(image, start, end, newName) {
// Compute the bits we need to extract.
var pattern = 0;
for (var i = start; i <= end; i++) {
pattern += Math.pow(2, i);
}
// Return a single band image of the extracted QA bits, giving the band
// a new name.
return image.select([0], [newName])
.bitwiseAnd(pattern)
.rightShift(start);
};
// A function to mask out cloudy pixels.
var cloud_shadows = function(image) {
// Select the QA band.
var QA = image.select(['BQA']);
// Get the internal_cloud_algorithm_flag bit.
return getQABits(QA, 7,8, 'Cloud_shadows').eq(1);
// Return an image masking out cloudy areas.
};
// A function to mask out cloudy pixels.
var clouds = function(image) {
// Select the QA band.
var QA = image.select(['BQA']);
// Get the internal_cloud_algorithm_flag bit.
return getQABits(QA, 4,4, 'Cloud').eq(0);
// Return an image masking out cloudy areas.
};
var maskClouds = function(image) {
var cs = cloud_shadows(image);
var c = clouds(image);
image = image.updateMask(cs);
return image.updateMask(c);
};
var l8_mayo_free = l8_mayo.map(maskClouds);
Map.addLayer(l8_mayo_free, visParams, 'l8_mayo collection without clouds');
Link: https://code.earthengine.google.com/44cb99f8210800462da56640a275e674
There are some areas with cloud borders in image collection. Is a very cloudy area. BQA isn't perfect
In Arcgis, there is a way to offset graphically the contour line of a polygon like below. I'm trying to do the same thing in a styling method (i.e the solution can't be by buffering physically the layer). QGIS does have an offset option, but it's an XY type offset that moves the position of the line.
ANSWERS :
Way 1: use the geometry generator fill style. Just get in the box buffer($geometry,[your buffer distance]). The only difference is that Arcgis offsets the line whereas this solution creates a buffer, thus rounding the edges...
Way 2 : exact replica of Arcgis, inspired below by @iant answer : create in the style another style "Outline : Simple line". In this style, there is the offset box i was looking for ! It does exactly the same as in arcgis.
Answer
Another way to do this is to duplicate the layer and use a negative offset to go outside the polygon when styling with a simple line fill. Or just have a simple fill and a simple line with the offset.
I have joined two polyline shapefiles together that represent two months worth of data. I'm now trying to find any changes to the shape of the records between the two months (added vertices, deleted vertices).
Is this possible to detect these changes with a definition query?
If not, is there an efficient way to detect changes in the geometry of the record?
I have a "LENGTH" field and I am able to query records that have a different length but that doesn't necessarily give me all of the records who have had their vertices changed.
Answer
My direct Answer to your Question is "No - but Spatial Definition Queries are the subject of an ArcGIS Idea that you may want to vote for".
However, I think a much better approach is to use Feature Compare.
Alternatively, if you are more interested in a database approach to tracking changes through delta (add/delete) tables you may want to investigate Geodatabase Versioning.
How does one break the execution of a long-running process in the python console of ArcMap?
For example, how do you break a long iteration like the one below once it is started?
for i in range(1, 50):
# some long process which takes over a minute to complete
I want to stop the process at any time I want by pressing a keyboard combination, clicking a button or something similar while the code I have input in the console window is running.
Answer
Just tested it and ESC doesn't work either. ArcGIS just freezes for a moment and then continues. There doesn't seem to be a way to do force quit it once it runs in the ArcGIS Python console. You can't kill it using Task Manager either as the Python process doesn't show up there.
If you really want to be able to force quit it, you might want to consider using an IDE like IDLE rather than doing it in the ArcGIS Python console. In IDLE for example, you can use Ctrl + Z or Ctrl + C to terminate the execution.
I am fairly new to python and seek guidance for a question which might sound trivial to many.
Is there a way to use 'wget' in a python script to download raster files from a server and process them in the same script?
Answer
Python has urllib2 built-in, which opens a file-pointer-like object from a IP resource (HTTP, HTTPS, FTP).
import urllib2, os
# See http://data.vancouver.ca/datacatalogue/2009facetsGridSID.htm
rast_url = 'ftp://webftp.vancouver.ca/opendata/2009sid/J01.zip'
infp = urllib2.urlopen(rast_url)
You can then transfer and write the bytes locally (i.e., download it):
# Open a new file for writing, same filename as source
rast_fname = os.path.basename(rast_url)
outfp = open(rast_fname, 'wb')
# Transfer data .. this can take a while ...
outfp.write(infp.read())
outfp.close()
print('Your file is at ' + os.path.join(os.getcwd(), rast_fname))
Now you can do whatever you want with the file.
I am doing a household survey in Kenya, and I need to ensure that my survey team visits every house within a given area. I've been going to google maps and right-clicking on roofs to get coordinates via the "what's here" option. This is really an inefficient use of time.
I'd prefer to have a high-resolution image that I could upload to google earth engine, and then perform a machine-learning algorithm to identify roofs. I'd want to get a text file of the lat/longs of all the houses, which I could then cluster by distance and give to my surveyors, both in map and in checklist format.
But the images in google maps satellite view are copyrighted and not available in earth engine. So I need images at high-enough resolution to identify corrugated iron and straw thatch roofs, with roughly ~3km^2 spatial extent.
What are some good sources for such data? Ideally that will give it freely or cheaply to graduate students and/or nonprofits?
I have a series of maps representing a 1000m diameter around a point. These maps contain 40 different types of polygons that represent a land use (ex. Corn crop, Pasture).
I want to select all polygons that are within the 1000m boundary and find out each of their areas.
I've attempted to create unions and delete the exterior polygons, but I'm only limited to performing a union on 2 features which would make that task very tedious. Is there a faster and simpler way of performing this task?
My goal is to create an attbiute table, that will have a unique identifer for each polygon and then a column with the SHAPE_AREA that I can export to excel.
I created the maps using ArcGIS, I have access to all licenses and other GIS programs. All polygons are shapefiles, including the boundary circles.
Answer
I've determined the best possible route for my question.
I create a new field in the attribute table for each shapefile called "LandUse". I designated a number for each LandUse (I.e. Corn = 3)
I used the 'Merge' tool on all the LandUse shapefiles. Now I can use the Clip tool to determine Land Uses within different sized boundaries including 1000m (Boundaries as the 'clip features').
I use the field calculator in the Attribute table of the newly clipped shapefile to determine area and the summarize option to create a table of areas for each boundary)
Thanks everyone for your answers, they were all helpful!
I have a 10m DEM covering part of a state, a 20m DEM covering the whole state, and a 90m DEM covering the whole country. How can I merge them together, getting the best quality available in each area?
The 10m and 20m are provided in ADF format, but I've converted them to .tif using ogr2ogr.
Using gdal_merge doesn't work, because the boundaries of each area are irregular. This is the result:
(The northeastern edge of the state is irregular, but there obviously shouldn't be the gap in the middle.)
UPDATE If it's useful, I also have a shapefile for each 10m DEM, marking its extent.
UPDATE 2
Ok, here's what the whole state looks like with the 20m DEM.
The process I'm going through so far is this:
for f in dtm20m dtm10m_nw dtm10m_e; do
export GDAL_CACHEMAX=1000
echo -n "Re-projecting: "
gdalwarp -s_srs EPSG:3111 -t_srs EPSG:3785 -r bilinear ./vmelev_${f}/${f} ${f}-3785.tif
echo -n "Generating hill shading: "
gdaldem hillshade -z 5 $f-3785.tif $f-3785-hs.tif
echo and overviews:
gdaladdo -r average $f-3785-hs.tif 2 4 8 16
done
I am working for the Town of Yadkinville, NC and would like to know how to create a layer of directional flow arrows for sewer lines to be able to turn them on and off at will?
I am trying to create an new unique id field in a feature class. I already have one field called "SITE_ID_FD", but it is historical. The format of the unique value in that field isn't what our current format is, so I am creating a new field with the new format.
Old Format = M001, M002, K003, K004, S005, M006, etc
New format = 12001, 12002, 12003, 12004, 12005, 12006, etc
I wrote the following script:
fc = r"Z:\test.gdb\testfc"
x = 12001
cursor = arcpy.UpdateCursor(fc)
for row in cursor:
row.setValue("SITE_ID", x)
cursor.updateRow(row)
x+= 1
This works fine, but it populates the new id field based on the default sorting of objectID. I need to sort 2 fields first and then populate the new id field based on that sorting (I want to sort by a field called 'SITE' and then by the old id field "SITE_ID_FD")
I tried running the script from the python window in ArcMap after manually sorting the 2 fields in hopes that the python would honor the sort, but it doesn't. I'm not sure how to do this in python. Can anyone suggest a method?
I use QGIS 2.14.2. I set up labels names of a line vector and set it up as curved - looks like in this screenshot
But when I want to print it in print composer it doesn't show up 
It only shows labels if they aren't curved. 
What can I do to fix it. Curved labels are looking really nice (for compared to parallel one.
Answer
The problem is with map scale. Too short line for those words.
I want to mask pixels based on the quality band information for image collection Landsat 5. this is a code but the image-collection is empty afterward. I would be thankful if anyone could help me.
var startDate = '1992-01-01'
var endDate = '1993-01-01'
var cloudCover = 30
var collection1992 = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
.filter(ee.Filter.eq('WRS_PATH', 170))
.filter(ee.Filter.eq('WRS_ROW', 32))
.filterDate(startDate,endDate)
.filterMetadata('CLOUD_COVER','less_than',cloudCover)
.sort('system:time_start')
// cloud mask 66,130
var collection1992_noclouds = collection1992.map(function(img) {
var mask1 = img.select('pixel_qa').eq(66).eq(130)
return img.UpdateMask(mask1)
})
Map.addLayer(collection1992_noclouds, {}, 'collection1992_noclouds')
I am trying to publish some ecw files with geoserver 2.7.2 on ubuntu 15.04 x64. I downloaded gdal plugin, gdal_data and native libraries for image-io 1.1.11 and set environment variables; GDAL_DATA, LD_LIBRARY_PATH.
I can see new store types when adding new one (ECW, EHdr, JP2ECW etc.). But when I try to save an ECW store with a valid ecw file which I tested with ArcMap, geoserver crashes with this error:
terminate called after throwing an instance of 'std::length_error' what(): basic_string::_S_create Aborted (core dumped)
Do you have any idea about this problem?
Extremely new to QGIS and was hoping to merge 2 layers and maintain an "overlay" keeping layer transparency and color rendering?
I seem to loose the top layers layer properties when merging and it reverts to it's original look.
Before "Merge": 
After "Merge": 
Using QGIS 3.2.2
File geodatabase (fgdb) includes numerous file geodatabase tables. As far as I know they exist as dbf files, but are within a Database.gdb.
In ArcCatalog, the file pathway resembles C:\Users\...\Database.gdb\Stats_AA.
How to read all of these dbf files into R (a statistical software)? What is the correct pathway to supply? The function used is read.dbf (in the foreign package).
Variants of
test<-read.dbf(file="C:/Users..Database.gdb/Stats_AA.dbf")
and
test<-read.dbf(file="C:/Users..Database/Stats_AA.dbf")
don't work. What's the correct "form" of the file name to be used, or, do I need to export all of the file geodatabase tables into some other form or location?
Answer
A simple solution is to use Table to dBase (multiple) to export your tables (Right click FGDB > Export > To dBase (multiple). You can also use this tool to export attribute tables contained within FGDB feature classes. Just drag and drop tables and or feature classes into the tool and specify an output folder. Of course, then you can loop through the folder containing the new dBase files using R.
I am working on a script that will take an unprojected raster, and using a UTM grid (shapefile), will obtain the correct SRS (I have a shapefile where each zone is attributed with its numerical SRS) and perform a warp/projection using GDAL on that raster. It is a quicker way than figuring out the right zone, and then manually doing the warp.
Procedurally, I am grabbing the raster corners and constructing WKT geometry (polygon); the objective is to do a spatial join with the correct feature in the UTM grid. Here's what I have so far:
import os
import ogr
from osgeo import gdal
def getRasterBox(raster):
src = gdal.Open(raster)
ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform()
lrx = ulx + (src.RasterXSize * xres)
lry = uly + (src.RasterYSize * yres)
return (uls, uly, lrx, lry)
def findUTMzone(raster_box, utmgrid):
# do the stuff
utmzone_path = r'C:\temp'
utmzone_file = 'utm_grid.shp'
utmzone = os.path.join(utmzone_path, utmzone_file)
utm_ds = ogr.Open(utmzone)
utm_layer = utm_ds.GetLayer()
rasterdir = r'C:\temp'
raster_file = 'n38.tif'
raster = os.path.join(rasterdir, raster_file)
rasterdetails = getRasterBox(raster)
rasterWKT = 'POLYGON(({} {}, {} {}, {} {}, {} {}, {} {}))'.format(rasterdetails[0], rasterdetails[1], rasterdetails[0], rasterdetails[3], rasterdetails[2], rasterdetails[3], rasterdetails[2], rasterdetails[1], rasterdetails[0], rasterdetails[1])
Using the QuickWKT plugin for QGIS, I can see that I'm getting a good WKT ring for the raster bounds. However, I'm unsure how to perform the spatial join in a way that grabs the feature attribute containing the correct SRS code. How do I do this part (under def findUTMzone(raster_box, utmgrid):)? I have seen a few arcpy examples, but am looking for something open-source.
My task is very straightforward: a user draws a polygon and queries all parcel polygons. The query string extends the URL limit. So I am posting XML to WFS. However, GeoServer always returned 0 feature, but I knew the query polygon overlapped several tax parcels. If I passed a single point, it worked fine. Did I miss anything for sending a polygon? I used the OpenLayers buildGeometry to create the polygon GML. I can't figure out what's wrong.
Please Help.
Here is my code:
var gml = new OpenLayers.Format.GML();
var polygonGML = gml.buildGeometry.polygon.apply(gml,[taArray[0].geometry])
console.log(polygonGML);
var plgGMLString = (new XMLSerializer()).serializeToString(polygonGML);
var geo_req_url = proxyHost + "http://localhost:8080/geoserver/wfs";
var query_xml = ' + ''
+ 'geom '
+ plgGMLString
+ ' ';
$.post(
geo_req_url,
query_xml,
function(data) {
geoResults = data.features;
remove_selected_features();
draw_sel_features();
},
"application/json"
);
Here is the printout of plgGMLString.
How does one add a field to an ArcGIS feature class with a boolean data type? That is an attribute where the allowed value is only one of a pair such as 1/0, Yes/No, On/Off, Present/Not-present, etc.
A boolean data type is not available in the Data Type pick-list.
Answer
Many DBMSes including Oracle do not have a Boolean column type, so this may explain why there is no Boolean field type in ArcGIS either.
@Emily's suggestion of using a coded value attribute domain is a good one and I believe this is the ESRI-recommended best practice.
The only other suggestion I would have is to use a CHAR field of length 1, as this way you can use the slightly more descriptive "Y" or "N", or "T" or "F", rather than 1's or 0's, and, according to this article, in Oracle 1-character CHAR columns are actually more efficient than 1-digit NUMBER columns.
I would like to calculate the sum of values in one column, based on values in another column. To illustrate my question I uploaded the picture below. The second column consists of zipcodes, the third total sales. I would like to know what the total sales per zipcode are? So in this case:
What is the sum of column three for column two's value 1023? Output then should be 4 (0+1+1+2+0). The sum for 1024 then should be 11,5 following the same logic.
As this is quite a large dataset I would like to have the QGIS field calculator (or any other tool) calculate a new field which states the total sales in the zipcode. So far I can't seem to get it working and was hoping someone here could help me.
As you can see I have some missing values, is this a problem in calculation? Perhaps good to mention is that it is a table.
I have a table in a Oracle SQL database in which I have a polygon geometry and two more columns of latitude and longitude. What I need is to show the points of lat and lng in a layer in geoserver.
I have done similar thing using PostGIS but in now with Oracle I face an error. First I have create a query which selects the two coloumns, transforms them to the desired SRID (EPSG:2100), REPLACE the '.' with the ',' and CAST them to float (cause they are characters).
When I execute the query in the Oracle SQL developer it seems that it works correctly. Also if I copy and paste teh coords they are shown in the correct position on the map.
But when I se my query in geoserver in order to create a sql view layer, although the geometry column is recognized (also the SRID), when I click on the "compute BBOX" option I get a geoserver error.
As I said I did something similar with postGIS and it was working fine. I also tried to use the same table (but its existing geometry; the polygon) and it works also great.
So I guess it must be related with my query but I can not guess what can cause this "unexpected geoserver error".
Thats my query:
SELECT ID, CO_NAME, ADJUSTMENT_FACTOR, PORT_LOCATION, VDSL_ENABLED,
SDO_CS.TRANSFORM(SDO_GEOMETRY(2100,4326,SDO_POINT_TYPE(
CAST(REPLACE(CENTRAL_OFFICES.LAT, '.', ',') AS FLOAT),
CAST(REPLACE(CENTRAL_OFFICES.LNG, '.', ',') AS FLOAT),NULL),NULL,NULL),2100) AS geom
FROM CENTRAL_OFFICES
Answer
Let's assume that your table is like this:
create table central_offices (
id number primary key,
lat varchar2(20),
lng varchar2(20)
);
insert into central_offices (id, lat, lng) values (1, '39.018483', '22.9983436');
commit;
Your strings are stored with a point as decimal separator. When you use this string in a statement that expects a number, Oracle will automatically do the proper type casting and conversion, i.e. parse the string into a number.
BUT: that parsing happens in the context of the locale you use in your session. Since 2100 is a Greek system, I assume you are in Greece, and so use a Greek locale where the decimal separator is a comma. If you try converting the above string to a number, you will get the following error:
SQL> select to_number('39.018483') from dual;
select to_number('39.018483') from dual
*
ERROR at line 1:
ORA-01722: μη αποδεκτός αριθμός
Replacing ',' with '.' on the fly is a valid solution. A simpler approach is just to override the locale for your session:
alter session set nls_numeric_characters = '.,';
select id,
sdo_cs.transform (
sdo_geometry (2001, 4326, sdo_point_type(lng, lat, null), null,null),
2100
).get_wkt()
from central_offices;
ID SDO_CS.TRANSFORM(SDO_GEOMETRY(2001,4326,SDO_POINT_TYPE(LNG,LAT,NULL),NULL,NULL)
---- -------------------------------------------------------------------------------
1 POINT (413136.397473566 4319017.56676377)
I'm using ArcGIS 10.3 (Desktop), and I'm trying to measure the distance between points following a path determined by a set of polylines. Hand measurement is not a viable option-- there are ~7500 points, and my polyline set is the global coastline. The ultimate goal is to get the spacing of these points along the polylines, so that I can analyze this with respect to a few other variables.
I can read python-- with a bit of effort-- but not write or edit it, so code-based solutions will not be particularly useful to me.
I have tried using the linear referencing tools-- I snapped the points to the lines, then transformed the polylines to routes and used 'locate features'-- but because my points are not (entirely) ordered along the path, using the difference between linear reference locations of subsequent points yields some fairly unlikely answers. In addition, spot checks of (hand-measured, using the built-in path length measuring tool) point-to-point along-polyline distance yield different answers than given by the corresponding linear references for the point pair. (Some of the answers were close enough that they could very well be measurement error on my part, but others were nearly 20%, despite my having zoomed in and snapped the measurement cursor to the polyline, proceeding vertex-to-vertex.)
Sorting first by linear distance seems to fix the first problem, but in light of the second one, I remain skeptical of this process.
I then tried simply splitting the polylines into segments ('split line at point') using the (snapped) points as the dividers, with the eventual goal of finding the segment length with the calculate geometry feature. This produced gibberish. I'm not sure what criteria ArcGIS used to split the lines, but it wasn't the point set I provided. Segments in the resulting shapefile ran right through tens of points, only to terminate at a place with no point. (Also, there were segments utterly empty of points...you get the gist.) I reiterate that the pointset was specifically snapped to the polylines, so this problem shouldn't be caused by excessive distance between the point and the polyline.
I'm at a bit of a loss, here: could someone point out where I might have gone astray and/or suggest a solution? Or, failing that, a new method that I could try? It seems like a fairly basic (not necessarily simple, but fundamental) analysis to do with geographic data, and I don't understand why it has presented so many difficulties.
And I don't care about format--worry about that later! I found some maps in Google maps, usually per state. I also have one map made my someone else, that no longer appears in the Google Maps search for user content. Maybe that user made it private, but I still can see it.
The reason I ask is related to my google maps navigation question on android.se. I think the easiest way to get this info on the droid is to create or use a shared map in Google Maps. I can find Florida, but no others show up.
Maybe it's a vocabulary issue? What else are are rest stops called?
Answer
I think you're going to have to hunt this down on a state-by-state basis. A good place to start is the state spatial data clearinghouses. Google "spatial data clearinghouse" and most states will show up. Unfortunately some states (like my home state) are years behind in things like GIS, so you may not find universal coverage.
Have you seen the mashup at Programmable Web? There's a contact page and you might be able to find out where they got the data they used. Doesn't look like the mashup works right now, though.
Another place to check is a page MIT has of spatial data sites at Online spatial data by state.
This won't help but I thought it was funny that someone did this: web site listing rest areas by highway and exit number.
If you are open to purchasing data, there's NAVTEQ and Tele Atlas, but I think you can probably access most of their data through some map service somewhere.
I haven't used a POI file before, but you can get one at POI directory.
I am trying to clip a .tif file using the extent of a polygon (.shp). I am using the conventional Raster/Extraction/Clipper tool in QGIS 2.18.10. After choosing the appropriate settings, I get the following message:
ERROR 1: Attempt to create ADRG dataset with an illegal data type (Float32), only Byte supported by the format.
I am asking the exact same question as uybfi in this thread, which was marked as a duplicate for no apparent reason. The referred answer fails to address the issue. My raster is already in a supported format (Float32), and when changing the pixel depth with gdal_translate to other supported formats (e.g. Int32), I still get the this error.
Any thoughts?
Answer
The issue had to do with QGIS wrongfully inputing -of ADRG instead of -of GTiff in the command line when choosing the output file destination. Below is a screenshot of the clipper window properly set-up to successfully perform the clip operation. Thanks to @MichaelStimon for resolving this.
I am trying to create a map of forest type diversity with ArcMap. I have a polygon forest layer with the different forest types as attributes (e.g. 0815, 0862, 0813, ..., which stands for 'mixed forest', 'oak forest', and so on.)
Using Patch Analyst I have created a Hexagon map and intersected it with the forest polygon layer. Usually I would now use 'field statistics' to calculate the variance of forest types for each Hexagon field and then join the result table to the Hexagon polygon layer to create a diversity map. Calculations don't work, however, since I have qualitative Features instead of quantitative ones.
Does somebody have an idea how to do this?
Answer
I solved the Problem:
--> this gives me what I want, in the new table there is a Count_ field that Displays not the number of Polygons, but the number of different forest types in the hexagon
I'm working on rasterizing some polygons using gdal_rasterize. This works great, but I want the rasterize command to place features with lower FIDs on top of those with higher FIDs. (Gdal_rasterize currently takes the latest polygon and burns that value onto the raster).
I can think of a way to do it by manually swapping the rings for the polygons, but is there a faster way reverse the order of features in a layer?
Hoping for a simple reverse sort function? I'm working with OGR/GDAL on Python 2.x
Answer
If you are happy using the ogr2ogr utility you can accomplish this feature reversal by including an SQL statement in the command. For example, the following will create a new shapefile with the features in reverse order:
ogr2ogr -sql "SELECT * FROM layer ORDER BY fid DESC" destfile.shp srcfile.shp
If you need to have this reversal happen inside of your Python script you can use the datasource's ExecuteSQL() method to return a new layer that can be passed to the gdal.RasterizeLayer() function:
rev_lyr = ds.ExecuteSQL("SELECT * FROM layer ORDER BY fid DESC")
gdal.RasterizeLayer(out_ds, [1], rev_lyr)
Does ArcGIS 10 have a Geoprocessor Programming Model by using Arcpy? is it possible to use ArcGIS 9.3 Geoprocessor Programming Model within ArcGIS 10?
I'm new to the leaflet/esri subject so please forgive me if it's just a simple question.
I wrote my first webmap-app. It uses different sources. The main aim is to show two basemaps hosted by MapBox with an AGS layer coming from a different source.
Mapbox does not make any trouble. The problem is the AGS layer. Its address is:
http://geoportal1.stadt-koeln.de/ArcGIS/rest/services/Stadtplanthemen/MapServer/9
My script shows four different behaviours. I tested it with two different AG Servers:
Test with http://geoportal1.stadt-koeln.de/ArcGIS/rest/services/Stadtplanthemen/MapServer/9 a: It shows up in IE10 correctly b: Not showing up in FF 33
Test with the same script but a different AGS http://services.arcgis.com/rOo16HdIMeOBI4Mb/arcgis/rest/services/Heritage_Trees_Portland/FeatureServer/0 a: Shows up in IE 10 b: Shows up in FF 33
Well, what is the error in this script? The geoportal1's AGS seems to work properly. FF understands the script. But what's not working? If I use L.esri.dynamicMapLayer instead of L.esri.featureLayer the layer shows up. But for further uses I need the L.esri.featureLayer.
I also published it here:
https://github.com/Esri/esri-leaflet/issues/412
Hopefully someone can help!
A simple map
I am trying to work with DEM rasters from the USGS using the raster calculator. I have tried to change the units from meters to feet using an expression such as:
"ned10m45111h8@1" * 3.28
However this returns NAN values for every cell. I have also tried it without the quotes around the raster name as shown in this website: http://spatialgalaxy.net/2012/01/25/using-the-qgis-raster-calculator/
Similarly, I have tried the mask code:
("ned10m45111h8@1" <= 1328.96)*"ned10m45111h8@1"
With and without quotes and get NAN returned for all values. Is there some default setting that I need to add, or some fundamental mistake I am making? I'm new to QGIS, but not to GIS and could easily do this with ARCGIS and spatial analyst.
From a given coordinate and radius I need to create a circle with 12 sectors. The sectors should start from north (0 degrees) in 30 degree steps clockwise. I could also imagine to create the sectors by triangles, if that would be easier.
Answer
I wrote a small script that generates polygonal sectors (for now it exports the data as GeoJSON, but you should easily be able to export the data as a Shapefile etc.). You can easily adjust the parameters as you can see on the screenshot:
from shapely.geometry import Point, Polygon
import math
import geojson
# initial parameters for segmentation
steps = 90 # subdivision of circle. The higher, the smoother it will be
sectors = 12.0 # number of sectors in the circle (12 means 30 degrees per sector)
radius = 90.0 # circle radius
start = 315.0 # start of circle in degrees
end = 45.0 # end of circle in degrees
center = Point(23,42)
# prepare parameters
if start > end:
start = start - 360
else:
pass
step_angle_width = (end-start) / steps
sector_width = (end-start) / sectors
steps_per_sector = int(math.ceil(steps / sectors))
# helper function to calculate point from relative polar coordinates (degrees)
def polar_point(origin_point, angle, distance):
return [origin_point.x + math.sin(math.radians(angle)) * distance, origin_point.y + math.cos(math.radians(angle)) * distance]
features = []
for x in xrange(0,int(sectors)):
segment_vertices = []
# first the center and first point
segment_vertices.append(polar_point(center, 0,0))
segment_vertices.append(polar_point(center, start + x*sector_width,radius))
# then the sector outline points
for z in xrange(1, steps_per_sector):
segment_vertices.append((polar_point(center, start + x * sector_width + z * step_angle_width,radius)))
# then again the center point to finish the polygon
segment_vertices.append(polar_point(center, start + x * sector_width+sector_width,radius))
segment_vertices.append(polar_point(center, 0,0))
# create feature
features.append(geojson.Feature(
geometry=Polygon(segment_vertices))
)
# prepare geojson feature collection
res = geojson.FeatureCollection(
features = features
)
# write to file
f = open('sector.json', 'w')
f.write(geojson.dumps(res))
f.close()
I've got fresh installation of ArcGIS 10.5 on my machine and was wondering from where and how can I get 64-bit background geoprocessing.
Help pages are not explaining anything regarding installation.
In previous versions it was (I think) located on the installation media (similarly as pointed to in this or this questions for 10.2):
My IT department informed me that this is not the case any longer and in some newer thread on geonet forum I found information that for version 10.3 this is downloadable from my.esri.com
(most likely with an institutional account..) but I cannot find any details regarding 10.5 version.
Has anyone managed to get and install it?
Answer
It's downloadable for 10.5 as well, but you need to have access to the Downloads page for your organization on my.esri.com. Most likely, someone in your organization can either download it for you or they can grant you the ability to download it for yourself. I'm pretty sure I remember this download, and most other addons/patches for ArcGIS Desktop being freely accessible in the past (without signing in to an Esri account).
Here it is on my downloads page. The filename is ArcGIS_Desktop_BackgroundGP_105_154031.exe.
I am running a brand new out of the box Samsung Galaxy Tab S2 with Android 6 and trying to get QGIS QField to operate correctly. So far I have created a QGIS project in QGIS desktop, packaged it with QfieldSync and uploaded it to the tablet, but I'm struggling to do anything with it.
It opens and I can see the individual layers in 'Digitise' but only the layer that I copied is visible in the data frame, the layers that I wish to edit and that are in 'offline editing' configuration are completely invisible. Is this normal (I am used to ArcGIS collector - where you can edit data and see it at the same time) ? I'd like to be able to see the existing data in these layers while I edit it. Is this possible and have a made an error or is there some other way to do this?
I am working with two polygon layers (counties & states) which I will 'intersect' to find which counties are within which states. Due to data inaccuracies I have to change all county polygons and make them a bit smaller. At this stage it doesn't really matter by how much, as long as I can do this change automatically.
I understand that polygons which are currently perfectly aligned now, will afterwards have a 'gap' between them, but that's fine.
I've tried 'Simpliyfy Geometries' but this only 'smoothens' edges and doesn't help in this case.
Is there any option in QGIS to reduce the size of polygons automatically? (i.e. move all polygon points 200 meters closer to the 'polygon mid-point').
I have a road network with 8.5 million edges but when I run pgr_createTopology it processes edges at a rate of about 1000 edges per second up until 360,000 edges where it slows down and eventually fails with the following error.
server closed the connection unexpectedly
This probably means the server terminated abnormally
before or while processing the request.
The connection to the server was lost. Attempting reset: Failed.
I am running postgres 10.1 postgis 2.4 pgrouting 2.5. What would like cause it to fail like this? Or is there any way that I could process the network in 100,000 edge chunks?
Answer
The following is what I am using. Some of it is specific to our deployment environment since we are using docker and some bash scripts to deploy and set up the server. You could easily get rid of all the argeparse/os.getenv and hardcode the connection if you wanted.
import argparse
from os import getenv
import psycopg2
parser = argparse.ArgumentParser()
parser.add_argument("-H", "--host", help="host location of postgres database", type=str)
parser.add_argument("-U", "--user", help="username to connect to the database", type=str)
parser.add_argument("-d", "--dbname", help="database name", type=str)
parser.add_argument("-p", "--port", help="port to connect to postgres", type=str)
args = parser.parse_args()
password = getenv('POSTGRES_PASSWORD')
conn = psycopg2.connect(
f"dbname={args.dbname} user={args.user} host={args.host} port={args.port} password={password}"
)
cur = conn.cursor()
print("connected to database")
cur.execute("SELECT MIN(id), MAX(id) FROM ways;")
min_id, max_id = cur.fetchone()
print(f"there are {max_id - min_id + 1} edges to be processed")
cur.close()
interval = 200000
for x in range(min_id, max_id+1, interval):
cur = conn.cursor()
cur.execute(
f"select pgr_createTopology('ways', 0.000001, 'the_geom', 'gid', rows_where:='id>={x} and id<{x+interval}');"
)
conn.commit()
x_max = x + interval - 1
if x_max > max_id:
x_max = max_id
print(f"edges {x} - {x_max} have be processed")
cur = conn.cursor()
cur.execute("""ALTER TABLE ways_vertices_pgr
ADD COLUMN IF NOT EXISTS lat float8,
ADD COLUMN IF NOT EXISTS lon float8;""")
cur.execute("""UPDATE ways_vertices_pgr
SET lat = ST_Y(the_geom),
lon = ST_X(the_geom);""")
conn.commit()
I hope this question is appropriate here and I apologize if it is more of a discussion question than a specific technical problem.
I have to buffer 150,000 polygons for a real estate project. While it's not impossibly slow, it got me wondering about 'the fastest way to geoprocess'. What packages, tips, or tricks could be used to make this faster? Does anyone have an idea how arcmap, arcpy, shapely, postgres compare in performance time?
EDIT: After asking around, I heard 'arcpy is slow!' a lot.
I have a point dataset which I'd like to Krige, ideally using an open-source software package. If possible, I'd also like to choose the semi-variogram model during the process to improve the estimation.
Answer
Depending on which Kriging type you want to apply, there are different packages to choose from:
The most common version is implemented for example in:
Simple Kriging uses the average of the entire data set while Ordinary Kriging uses a local average. Therefore, Simple Kriging can be less accurate, but it generally produces "smoother" results. It's implemented in:
Universal Kriging allows for consideration of drift in data. Implementations are included in:
Other Kriging Types
GRASS v.krige also supports Block Kriging.
HPGL implements a big number of less known Kriging methods (check the manual for more information on those):
SAGA offers different versions of both Ordinary and Universal Kriging.
Gstat krige additionally supports Block and Point Kriging.
I have a raster of continuous numeric values. A certain portion of said raster overlaps with a polygon of interest. I want to know which cells overlap, and assign values to the raster's attribute table corresponding to whether or not there is overlap (e.g. 0 for no overlap with polygon, 1 for overlap).
Ideally it would be great to be able to set an overlap threshold so that only raster cells where there is more than 50% overlap are assigned a "1". Not a must, but would be very helpful.
I'm currently trying to create polygons that represent the areas surrounded by roads. The roads are stored in PostGIS as points / lines so I have flexibility over what I can do with them.
I'm basically trying to turn the white areas of this example line output into polygons:
Any ideas?
Here is the PostGIS to achieve this (assuming you have a road table full of lines):
SELECT (ST_Dump(ST_Polygonize(roads.geom))).geom AS the_geom FROM
(SELECT ST_Transform(ST_SetSRID(geom, 27700), 4326) AS geom FROM road_lines) AS roads
Answer
Some hints:
I am having trouble exporting MXDs to PDF. If I export and MXD to PDF some point symbols appear as black dots (not the colored dots they should be) and others do not. This happens in different data sets in the same map.
I have embed fonts and convert marker symbols checked in the Export to PDF tool. I do not use non-USA thousand separators.
I have read this ESRI literature about diagnosing export problems:
http://support.esri.com/es/knowledgebase/techarticles/detail/17783
I had the IT department rebuild my Citrix profile, reinstall the PDF drivers, and I did a reinstallation of ArcGIS 10.3 for Desktop. The first map exported after all this work looked fine, but all subsequent exports revert back to the black dot problem.
I have point of data in mysql server. I am displaying this layer as WMS using OpenLayers and GeoServer. How can I add label to this layer?
Answer
For labeling in WMS from geoserver, you will have to define
circle
#FF0000
6
![]()
#000000
For more complex examples go through this link:
http://docs.geoserver.org/stable/en/user/styling/sld-cookbook/points.html#point-with-styled-label
I got the following error message after trying to import cloud-free Sentinel-2 imagery for a polygon area covering the entire Tiwi Islands, Australia.
true-colour image: Layer error: internal error
// Filter to only include images intersecting Tiwi.
var polygon = ee.Geometry.Polygon({
coords: [[
[129.90415998438357,-12.027354108834277], [131.63999982813357,-12.027354108834277],
[131.63999982813357,-11.101721188177093], [129.90415998438357,-11.101721188177093],
[129.90415998438357,-12.027354108834277]]],
geodesic: false
});
// Define the image collection we are working with by writing this command
var image = ee.Image(sent2
// filter to get only images in the date range we are interested in
.filterDate("2018-09-01", "2018-11-30")
// sort the collection by a metadata property, in our case cloud cover is a very useful one
.sort("CLOUD_COVERAGE_ASSESSMENT")
// select the first image out of this collection - i.e. the most cloud free image in the date range
.first());
// print the image to the console.
print("A Sentinel-2 scene:", image);
// Define visualization parameters in a JavaScript dictionary for true colour rendering. Bands 4,3 and 2 needed for RGB.
var trueColour = {
bands: ["B4", "B3", "B2"],
min: 0,
max: 3000
};
// Add the image to the map, using the visualization parameters.
Map.addLayer(image, trueColour, "true-colour image");
What I've got:
1: A vector-layer with administrative areas of my country
2: A raster-layer with the slope of the terrain of the same country (it is derived from SRTM elevation data with QGIS)
What I want:
A new vector-layer with the 'flat' land of all administrative areas.
How do I achieve this in QGIS?
A) How to extract the area below a certain slope (in the raster-layer)? - The option for rule-based display is not available (maybe this is just there for vector-layers).
B) How to 'intersect' the flat-land-raster-layer(-mask) with my admin-vector-layer?
Update:
A) Actually the "Clip to MinMax" with an appropriate set Min- and Max-value (e.g. Min: 10, Max: 100 - or inverse Min: 0, Max: 10.1) does the job. Now I want to use this as a 'clipping-mask' for my vector-layer, but how?
Answer
In QGIS you can use the gdal_polygonize.py function available through the menu Raster > Conversion > Polygonize (raster to vector). Once you have created the new polygon shapefile, then filter the data by slope. Next, delete the non-"flat" features from the shapefile by selecting them based on the slope values.
You can use the Vector > Geoprocessing Tools > Symmetrical difference tool to cut the remaining slope polygons along political boundary lines.
Of course, there are likely to be places where steep terrain is present in an otherwise predominantly flat region. I think your definition of "vector layer with the 'flat' land of all administrative areas" covers this issue, though.
There is a handy video on youtube about converting a raster to a set of polygons.
I need to perform a simple calculation in ArcGIS 10.1 Raster Calculator. I am using a raster layer that contains many cells that have a value of zero. I need to exclude these cells in the expression because they are skewing the results.
Here is the simple raster calculator expression that I am using: "rho"* 1005*(((("raster")*.02)-272.15)-"temp")/30)
I need to add to this expression so that only values greater than zero in the raster layer "raster" are consider in the expression.
Answer
You could try setting your Zero values to NULL or "NODATA". Then as I have understood it, NODATA values will not get processed. You can use the SetNull tool.
I have a huge amount of data in the form of CSV containing UTM coordinates as X and Y and an elevation value as Z information. I need to convert these data into a DEM as GeoTIFF for further analysis. In this case, a huge amount means 16 m. lines, with one point in X, Y and Z per line. The points are equally distributed, therefore an interpolation is not needed; each point just needs to be converted into a raster cell.
The original data came without separator, with fixed column widths. I already figured out how to convert the file syntax to use a separator instead of fixed widths and eliminate all space characters, using the stream text editor sed. From here on, normally my workflow would be to import the data into ArcGIS by creating a feature class from the X, Y and Z data and as a second step, converting the point shapefile into a GeoTIFF, using the Point to Raster tool. However, the file I currently have is way too big for this process.
Instead of the above described workflow, I was looking for an efficient alternative and discovered GDAL. However, in gdal_translate, the closest supported format I can find in the supported filetype list, is ASCII grid but no comma separated XYZ. Another difficulty is, that I have UTM coordinates, while most examples seem to use decimal degree coordinates. However, I need to stay within the UTM system (or at least, my output GeoTIFF needs to be in a UTM coordinate system).
So I am looking for a way to convert the CSV XYZ into a GeoTIFF, using GDAL, but so far wasn't able to find examples dealing with this exact problem. I would be very happy for some hints or even code examples.
Answer
You can do this using GDAL, it directly supports XYZ format. It doesn't matter if your coordinates are UTM, gdal_translate will output in the same coordinate system.
So to convert to GeoTIFF is as simple as:
gdal_translate test.xyz test.tif
Look at the GeoTIFF doc for output options (such as compression) and the gdal_translate doc for more usage info. In particular, you should specify what the coordinate system is with the -a_srs parameter.
-a_srs srs_def:
Override the projection for the output file. The srs_def may be any of the usual GDAL/OGR forms, complete WKT, PROJ.4, EPSG:n or a file containing the WKT.
gdal_translate -a_srs EPSG:12345 test.xyz test.tif
Comma/space separated and fixed column widths, with and without a header row are supported.
The supported column separators are space, comma, semicolon and tabulations.
$ head -n 2 test_space.xyz
x y z
146.360047076550984 -39.0631214488636616 0.627969205379486084
$ gdalinfo test_space.xyz
Driver: XYZ/ASCII Gridded XYZ
Files: test_space.xyz
Size is 84, 66
Coordinate System is `'
Origin = (146.359922066953317,-39.062997159090934)
Pixel Size = (0.000250019195332,-0.000248579545455)
Corner Coordinates:
Upper Left ( 146.3599221, -39.0629972)
Lower Left ( 146.3599221, -39.0794034)
Upper Right ( 146.3809237, -39.0629972)
Lower Right ( 146.3809237, -39.0794034)
Center ( 146.3704229, -39.0712003)
Band 1 Block=84x1 Type=Float32, ColorInterp=Undefined
Min=0.336 Max=0.721
$ head -n 2 test_commas.xyz
x, y, z
146.360047076550984, -39.0631214488636616, 0.627969205379486084
$ gdalinfo test_commas.xyz
Driver: XYZ/ASCII Gridded XYZ
etc...
$ head -n 2 test_formatted.xyz
x y z
146.3600471 -39.06312145 0.627969205
$ gdalinfo test_formatted.xyz
Driver: XYZ/ASCII Gridded XYZ
etc...
The only gotchas I'm aware of are:
The file has to be sorted correctly (by Y, then X).
Cells with same Y coordinates must be placed on consecutive lines. For a same Y coordinate value, the lines in the dataset must be organized by increasing X values. The value of the Y coordinate can increase or decrease however.
$ head -n 5 test.csv
x,y,z
146.3707979,-39.07778764,0.491866767
146.3787985,-39.07157315,0.614820838
146.3637974,-39.07132457,0.555555582
146.3630473,-39.07579901,0.481217861
$ gdalinfo test.csv
ERROR 1: Ungridded dataset: At line 3, too many stepY values
gdalinfo failed - unable to open 'test.csv'.
$ tail -n +2 test.csv| sort -n -t ',' -k2 -k1 > test_sorted.xyz
$ head -n 5 test_sorted.xyz
146.3600471,-39.07927912,0.606096148
146.3602971,-39.07927912,0.603663027
146.3605471,-39.07927912,0.603663027
146.3607971,-39.07927912,0.589507282
146.3610472,-39.07927912,0.581049323
$ gdalinfo test_sorted.xyz
Driver: XYZ/ASCII Gridded XYZ
etc...
In PyQGIS, how can I get the number of features that have been created in a given layer and that are pending to be saved in the layer when "Toggle editing" or "Saving layer edits"?
Answer
Say that you have a vector layer referenced:
lyr = iface.activeLayer()
At this point, I assume you start an edit session and digitize some features.
Now you can use the QgsVectorLayerEditBuffer class, in this way:
if lyr.editBuffer():
print len( lyr.editBuffer().addedFeatures() ), "features to add!"
addedFeatures gives you a dictionary of new features that are pending. As you're interested in the count, you can use len().
In OpenLayers, is there an event which fires when the map starts drawing any layer, and another for when the map has finished drawing all layers?
I wish to display a "waiting" animation while any layer is still drawing.
The help file for LoadingPanel mentions something which sounds promising, but the sample doesn't work and it mentions tiles (I'd like it to work on vector layers too).
I'm looking for the equivalent of dojo.connect(map,"onUpdateEnd"... in the ArcGIS JS API:
Fires after layers that are updating their content have completed. This event is often used in combination with onUpdateStart to display a "Map is busy" or "Loading? " message as visual feedback to the end-user.
(As a workaround, I could count the layers in the map, and monitor the status of each layer using layer.events.register("loadend"... until they have all finished. Are there other options?)
Thanks
Answer
A loading panel is available as an add-in and was first mentioned here. You can either use that or hack it to suit your needs; anyway imho it provides good enough functionality.
The linked example does not work so I took the liberty of assembling a working jsfiddle:
Are there any tools that can take a DTM, in asc or geotiff format and generate a mesh from this that can used in meshlab or similar software.
I have asc DTM file that I need to convert to a ply mesh file. Is this possible?
In response to the comment. Then these DTM asc files are not the ones that meshlab understands.
The files starts with:
ncols 6250
nrows 6250
xllcorner 630000.000000000000
yllcorner 6070000.000000000000
cellsize 1.600000000000
NODATA_value -9999
Answer
A DTM raster can be represented by triangle meshes by finding a set of non-overlapping triangles that covers the entire mesh and approximates the elevation field. There are two different types of triangle meshes that can be used for this purpose:
In most applications, if you need to deal with elevation meshes, you'd go with a TIN since throwing out redundant or nearly redundant information allows for more efficient computations. However, creating TINs from rasters isn't straightforward, since there are many different triangulations that approximate a grid with the same error, but using different vertex sets.
Software for creating TINs
Raster to TIN function from the 3D Analyst toolbox.Grid to TIN (Surface Specific Points) function, followed by Export TIN to Stereo Lithography File (STL) function to export the TIN to a mesh format readable by Meshlab.Software for creating TRNs
Grid to TIN function.Roll your own solution: Implementing it isn't particularly hard. Here's a Python script that uses the GDAL library to read a raster DTM, and then writes out a binary PLY mesh.
Save the script as gdal_rastertotrn.py, then call it using python gdal_rastertotrn.py .
Here's an example. Converting a raster DTM of Crater Lake called crater_lake.tif...
... by calling python gdal_rastertotrn.py crater_lake.tif crater_lake.ply, and opening the resulting crater_lake.ply in Meshlab:
Here's the (unpolished) script. Since it uses the GDAL library, it can convert all raster types supported by GDAL. It only writes PLY files.
#!/usr/bin/python
import sys
import numpy as np
from osgeo import gdal
def write_ply(filename, coordinates, triangles, binary=True):
template = "ply\n"
if binary:
template += "format binary_" + sys.byteorder + "_endian 1.0\n"
else:
template += "format ascii 1.0\n"
template += """element vertex {nvertices:n}
property float x
property float y
property float z
element face {nfaces:n}
property list int int vertex_index
end_header
"""
context = {
"nvertices": len(coordinates),
"nfaces": len(triangles)
}
if binary:
with open(filename,'wb') as outfile:
outfile.write(template.format(**context))
coordinates = np.array(coordinates, dtype="float32")
coordinates.tofile(outfile)
triangles = np.hstack((np.ones([len(triangles),1], dtype="int") * 3,
triangles))
triangles = np.array(triangles, dtype="int32")
triangles.tofile(outfile)
else:
with open(filename,'w') as outfile:
outfile.write(template.format(**context))
np.savetxt(outfile, coordinates, fmt="%.3f")
np.savetxt(outfile, triangles, fmt="3 %i %i %i")
def readraster(filename):
raster = gdal.Open(filename)
return raster
def createvertexarray(raster):
transform = raster.GetGeoTransform()
width = raster.RasterXSize
height = raster.RasterYSize
x = np.arange(0, width) * transform[1] + transform[0]
y = np.arange(0, height) * transform[5] + transform[3]
xx, yy = np.meshgrid(x, y)
zz = raster.ReadAsArray()
vertices = np.vstack((xx,yy,zz)).reshape([3, -1]).transpose()
return vertices
def createindexarray(raster):
width = raster.RasterXSize
height = raster.RasterYSize
ai = np.arange(0, width - 1)
aj = np.arange(0, height - 1)
aii, ajj = np.meshgrid(ai, aj)
a = aii + ajj * width
a = a.flatten()
tria = np.vstack((a, a + width, a + width + 1, a, a + width + 1, a + 1))
tria = np.transpose(tria).reshape([-1, 3])
return tria
def main(argv):
inputfile = argv[0]
outputfile = argv[1]
raster = readraster(inputfile)
vertices = createvertexarray(raster)
triangles = createindexarray(raster)
write_ply(outputfile, vertices, triangles, binary=True)
if __name__ == "__main__":
main(sys.argv[1:])
I am using ArcGIS 9.3 and python 2.5. My code is shown at the bottom of this post.
I have a raster containing fire intensity information at 200m cell resolution. I run the script and the first count cell part returns the size of the fire. This works fine.
I then clip the raster using a polygon representing a forested catchment. The new raster represents the area of the catchment that is burnt. When I try to use the same code to count the number of cells of the new raster, it always gives me an answer of zero, even if there is definitely a raster created by the clip tool.
I have tried every combination of every tool I can think of.
my code...
fire = str(ignitionID) + "_" + str(weatherID) + "_Intensity.asc"
print "fire: " + fire
# ---------------------------------------------------------------------------
# this section calculates the size of the fire for model calibration purposes
# ---------------------------------------------------------------------------
# Process: Build Raster Attribute Table...
the_fire = input_fires + fire
shutil.copy(the_fire, "C:\\FIRESTORM\\Workspace") # copy the fire to workspace
the_fire = "C:\\FIRESTORM\\Workspace\\" + fire
try:
gp.BuildRasterAttributeTable_management(the_fire, "NONE")
resultf = gp.GetCount_management(the_fire) # count num of rows in raster
countf = int(resultf.GetOutput(0))
fire_area = countf * 4 # each cell has an area of 4 ha
print "the fire area = " + str(fire_area) + " ha"
except: # the attribute table cannot be created coz the fire is too large
fire_area = 300000 # default value for large fires 300000 ha
print "the fire area EXCEEDS 262,140 ha"
# -----------------------------------------------------------------------------
# this section calculates the area of the catchments that are burnt by the fire
# -----------------------------------------------------------------------------
# Local variables...
catch_extent = input_catch + "WS_Catchment_region.shp" # catchment boundary
fire_raster = "C:\\FIRESTORM\\Workspace\\fire_raster"
fire_ascii = the_fire
Output_raster = "C:\\FIRESTORM\\Workspace\\Extract_ASCI1"
Output_ASCII = "C:\\FIRESTORM\\Workspace\\rastert_extract1.asc"
catch_burn_area = 0 # reset catchment burn area
# Process: ASCII to Raster...
gp.ASCIIToRaster_conversion(fire_ascii, fire_raster, "INTEGER")
# Process: Define Projection... based on projection of catch_extent
desc = gp.Describe(catch_extent)
spatialRef = desc.SpatialReference
gp.DefineProjection_management(fire_raster, spatialRef)
# Process: Extract by Mask...
gp.ExtractByMask_sa(fire_raster, catch_extent, Output_raster)
# Convert raster to ascii
gp.RasterToASCII_conversion(Output_raster, Output_ASCII)
# Process: Build Raster Attribute Table...
gp.BuildRasterAttributeTable_management(Output_ASCII, "OVERWRITE")
# count number of rows in raster
resulti = gp.GetCount_management(Output_ASCII)
counti = int(resulti.GetOutput(0))
catch_burn_area = counti * 4 # each cell has an area of 4 ha
print "catchment area burnt = " + str(catch_burn_area) + " ha"
I am looking to calculate the minimum distance from the Center of a polygon (A crescent shaped polygon would have a center inside the polygon) to its edge. I have parcel data and I am trying to calculate in my model I created how tall I can build a cell tower on that parcel. The maximum height that I can build the tower is equal to the minimum distance from the center of the parcel to its edge. Is their a simple way to calculate this? Once I have the minimum distance calculated I can complete my query for my model.
Answer
The following approach requires you to have an Advance license:
These two steps can be wrapped up in a model.
In a recent article by Stewart and colleagues in IJHG I stumbled upon interesting technique of visualizing data using ringmaps. Some more information referenced in the article about this technique here and here.
[Source]
From the article:
Ring maps were created in Adobe Illustrator through the application of a custom script that dynamically drew, distributed, and symbolized all graphic map elements. The values for symbolization were read from a Comma Separated Value (CSV) file that contained all county attribute data. Three county-level ring maps were developed.
Is there a way to automate implementation of a map like that in ArcGIS Desktop 10?
Answer
Recent article in ArcUser offers scripts that seem to be the closest thing. Although it doesn't use box plots, the code might be a good starting point to implement that! (via Matt Artz)
