texts:
New QGIS user, 2.6 Brighton.  Periodically my WMS links to USGS terrain (typical contour) and ortho-imagery just go bad and give me the following view.  This has happened before, but not for a while.  Any thoughts on this?  Can't get them to renew - just get the checkerboard image for both.  I'm online, all other web / other links working just fine. I've not done anything or altered the project in any form on last or current opening.

(24 hrs later...still no resolution to this...ortho-imagery link works just fine. USGS terrain will not show up.)


texts: How to get raster corner coordinates using Python GDAL bindings?

texts:
Is there a way to get the corner coordinates (in degrees lat/long) from a raster file using gdal's Python bindings?

A few searches online have convinced me that there is not, so I have developed a work around by parsing the gdalinfo output, it's somewhat basic but I thought it might save some time for people who might be less comfortable with python. It also only works if gdalinfo contains the geographic coordinates along with the corner coordinates, which I don't believe is always the case.

Here's my workaround, does anyone have any better solutions?

gdalinfo on an appropriate raster results in something like this midway through the output:

Corner Coordinates:
Upper Left  (  -18449.521, -256913.934) (137d 7'21.93"E,  4d20'3.46"S)
Lower Left  (  -18449.521, -345509.683) (137d 7'19.32"E,  5d49'44.25"S)
Upper Right (   18407.241, -256913.934) (137d44'46.82"E,  4d20'3.46"S)
Lower Right (   18407.241, -345509.683) (137d44'49.42"E,  5d49'44.25"S)
Center      (     -21.140, -301211.809) (137d26'4.37"E,  5d 4'53.85"S)

This code will work on files who's gdalinfo look like that. I believe sometimes the coordinates will be in degrees and decimals, rather than degrees, minutes and seconds; it ought to be trivial to adjust the code for that situation.

import numpy as np
import subprocess

def GetCornerCoordinates(FileName):
    GdalInfo = subprocess.check_output('gdalinfo {}'.format(FileName), shell=True)
    GdalInfo = GdalInfo.split('/n') # Creates a line by line list.
    CornerLats, CornerLons = np.zeros(5), np.zeros(5)
    GotUL, GotUR, GotLL, GotLR, GotC = False, False, False, False, False
    for line in GdalInfo:
        if line[:10] == 'Upper Left':
            CornerLats[0], CornerLons[0] = GetLatLon(line)
            GotUL = True
        if line[:10] == 'Lower Left':
            CornerLats[1], CornerLons[1] = GetLatLon(line)
            GotLL = True
        if line[:11] == 'Upper Right':
            CornerLats[2], CornerLons[2] = GetLatLon(line)
            GotUR = True
        if line[:11] == 'Lower Right':
            CornerLats[3], CornerLons[3] = GetLatLon(line)
            GotLR = True
        if line[:6] == 'Center':
            CornerLats[4], CornerLons[4] = GetLatLon(line)
            GotC = True
        if GotUL and GotUR and GotLL and GotLR and GotC:
            break
    return CornerLats, CornerLons

def GetLatLon(line):
    coords = line.split(') (')[1]
    coords = coords[:-1]
    LonStr, LatStr = coords.split(',')
    # Longitude
    LonStr = LonStr.split('d')    # Get the degrees, and the rest
    LonD = int(LonStr[0])
    LonStr = LonStr[1].split('\'')# Get the arc-m, and the rest
    LonM = int(LonStr[0])
    LonStr = LonStr[1].split('"') # Get the arc-s, and the rest
    LonS = float(LonStr[0])
    Lon = LonD + LonM/60. + LonS/3600.
    if LonStr[1] in ['W', 'w']:
        Lon = -1*Lon
    # Same for Latitude
    LatStr = LatStr.split('d')
    LatD = int(LatStr[0])
    LatStr = LatStr[1].split('\'')
    LatM = int(LatStr[0])
    LatStr = LatStr[1].split('"')
    LatS = float(LatStr[0])
    Lat = LatD + LatM/60. + LatS/3600.
    if LatStr[1] in ['S', 's']:
        Lat = -1*Lat
    return Lat, Lon

FileName = Image.cub
# Mine's an ISIS3 cube file.
CornerLats, CornerLons = GetCornerCoordinates(FileName)
# UpperLeft, LowerLeft, UpperRight, LowerRight, Centre
print CornerLats
print CornerLons

And that gives me:

[-4.33429444 -5.82895833 -4.33429444 -5.82895833 -5.081625  ]
[ 137.12275833  137.12203333  137.74633889  137.74706111  137.43454722]


texts: Static Tables and/or Tabs in ArcGIS Pro Layout

texts:
Is it possible to create a non-dynamic/static table in a layout in ArcGIS Pro?

I would like to create a table to act as a revision box. If not a table, can text be entered with tabs and tab stops instead?

The only help I can find about tables relates to dynamic tables that populate automatically based upon the map.


texts: VBScript, Error 999999 Syntax Error Calculate Field in ArcGIS

texts:
I am using a field [all scores] as a parameter agricultural classification. I tried using field calculator (VB script) with many kinds script, but it never works and giving me notification 99999, syntax error.

I use this pre-logic.

Dim output
    If [all_score]>45 and [all_score]<80 then
        output="Dryland"
    ElseIf[all_score]>79 and [all_score]<87 then
        output="Agricultural Wetland"
    ElseIf[all_score]>86 and [all_score]<137 then
        output="Plantation"
    Else
        output="Non- Agricultural"
    End If

Can someone tell me what I'm missing here?


texts: How to change field name in raster file?

texts:
I´m having a bit of an understanding problem with a DEM TIF file.  I have a DEM TIF file which is served with Geoserver.  When I click on the layer the WMSGetFeatureInfo gives me the following output:

GREY_INDEX
110.0

I realise that the 110.0 is the height, but what I don´t understand is why the raster field is called grey_index. The DGM was produced from a point shapefile which was then interpolated using spatialanalyst. The height field in the shapefile was called "z" so why is this then converted to "grey_index"?

For the WMSGetFeatureInfo I would like to change the field name to "Height" but how do I do this? Is this really a "Field" in the same sense as a shapefile dbf? Am I able to just change field names in raster files?

I have GDAL and ESRI libraries available if anyone can help.


texts:
I'm new here and was wondering if anyone could help?
I'm looking to make some 3D maps using DTM data but need to convert it into .STL format for my MasterCAM to read so that I can generate the CAD files needed.
Is this possible?


texts: Creating a link between QGIS and a spreadsheet

texts:
For a project, I am collecting data continuously, and will continue the datacollection for quite some time. I enter my data into an excel spreadsheet, and import it to QGIS 2.16 for display and further analysis (using geographical coordinates).
In order to be able to viasualize my results continuously, I have to regularily import my data to QGIS. It is timeconsuming, thus I am seaching for a way to make the process more effective.
Is there a way to connect my datasheet to QGIS, so that every time I update my datasheet, my QGIS map gets updates? Or does anyone know of a program that enables this?


texts: Merging and splitting raster files in QGIS

texts:
I have four .HGT topography tiles (tile1,tile2,tile3,tile4 in the photo below). They are in geographic coordinate system.
I would like to merge those four tiles, and then from that merged raster layer to cut off a piece (green rectangle in the photo) and to make a new raster layer from that cut off green rectangle piece.

The corners of the rectangle are suppose to be defined by adding 10 degrees of latitude and longitude to the central point (55,-2) of the green rectangle piece:

How can this be done in QGIS?
How can four .HGT raster files be merged and then a new raster layer cut off from it?


texts: ArcGIS Server 9.3.1 - Creating a MXD and Exporting a PDF From It

texts:
I've been tasked with developing a geoprocessing service in ArcGIS Server 9.3.1 that in addition to providing the resultant feature sets, it needs to produce a pdf in which the layers can be toggled on & off.  Upgrading to ArcGIS Server 10.0 is not going to happen (even with my begging & pleading).

The most simplistic workflow possible is outlined below.  Step 1 is easy, I am already doing this with the python geoprocessor - steps 2 and 3 are where it gets difficult.  I know creating a mxd and exporting to pdf are NOT possible with the python gp, so what are my options?  I have been looking at some SOE examples, but they seem far to complex for my needs.

Geoprocessing service does some magical processing and creates a geodatabase with result feature classes.
Create a mxd document containing result feature classes and apply symbology from layer files
Export pdf from mxd - store on server somewhere.

Just to clarify, the mxd is only temporary to create the pdf.  No services will need to come from it and it can be destroyed after pdf creation.  This workflow be initiated from a Flex app.


texts:
I am trying to establish a connection from ArcMap 10.0 -- Service Pack 5 to an Oracle DB (I thought it was 10g if might be 11g) on windows 7.. It is windows  have connected to the database via Catalog. However, I keep getting an error when trying to add query layer. I have attached a picture for more context. Any ideas as to why?
If I try to connect via oracle 11g it says failed to connect to the specified server. Underlying DBMS error.


texts: Projecting Geoseries Geopandas doesn't work


texts: Can any one tell me from where I can buy commercial High resolution Satellite images for Biodiversity research ? I need some websites that sells these kind of images ?


texts: Plot coordinates in ArcGIS 9.3 don't make sense

texts:
I have a layer of plot features in ArcGIS 9.3 that I want the UTM coordinates for. So I used the Feature to Point tool and then added XY coordinates to each point. I then viewed the Attribute table and noticed the coordinates don't make sense.

Here is the information about my map/projection:

Projected Coordinate System:    NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet
Projection: Lambert_Conformal_Conic
False_Easting:  2000000.00261667
False_Northing: 0.00000000
Central_Meridian:   -79.00000000
Standard_Parallel_1:    34.33333333
Standard_Parallel_2:    36.16666667
Latitude_Of_Origin: 33.75000000
Linear Unit:    Foot_US

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

I realize the units are in feet (and I ultimately want meters for UTM), so I converted my X and Y to meters by multiplying each value by 0.3048 in an excel spreadsheet.

The resulting units end up in the Pacific Ocean (they're supposed to be in NC in the U.S.).

Examples of my coordinates (from ArcGIS, so they're in feet):

X Coord              Y Coord
2008821.00994000000  823231.18663500000
2008796.70771000000  823457.50968900000

...converted to meters:

X Coord      Y Coord
612288.6438  250920.8657
612281.2365  250989.849

What am I missing here?

I'm an ArcGIS novice.


texts: ArcGIS 10.1 Network Analyst: Snapping stops from point-shapefile to route analysis layer doesn't work properly

texts:
I'm struggling with a task, which should be not that difficult....but anyway I do not find a solution.

Background
I'm map-matching tracks, generated by a GPS-device, to a road network. Therefore I make use of the ArcGIS 10.1 Network Analyst. I connected the GPS-Point measurements to polylines (1 for each track), created a buffer around them and made an intersection of the buffer-boundary with the street-network to determine point-barriers.

According to the paper of Ron Dalumpines and Darren M. Scott form McMaster University, Ontario, Canada I'd like to find the shortest path within the buffered area (confined by the generated point barriers).

Problem
The problem I'm struggling with is the following:
I'd like to define the stops (start- and end-point of the track) through a shapefile which contains the start- and endpoint for a certain track. The position of these points already coincidences with network-junctions points. Nevertheless ArcGIS seems not to be able, to link these points with the network and these "stop-points" get marked as "Unlocated".

If I use the default search tolerance of 5000 meters, ArcGIS is able to locate the stop-points....but the found ones are not at all the nearest ones to the position of the loaded shapefile.

Has anybody an idea what could cause this strange way of snapping? Has anybody already had a similar problem?

As a newbie I'm not allowed to upload pictures...otherwise I would have uploaded two pictures explaining my situation.

Thanks already in advance for your hints and advices ;)

Best Regards,

Dom

ETH Zurich, Switzerland


texts:
When I try to measure an element I get really big numbers,
I tried to measure by using in the attribute table -> open field calculator -> output field type (whole number (integer)) ->row_number - Geometry - area$


texts: Converting .kml to different projection in QGIS

texts:
I have run into a slight problem it would seem. Yesterday I downloaded qgis in the hope that it could convert a .kml overlay of some administrative boundaries from EPSG:3857 (Google maps) to EPSG:27700 (Ordnance Survey).

What I did: I loaded the .kml file into qgis as a vector file and then right click, save as and type in the new Projection I would like as in this photo:

However, when I load the newly created .kml in google earth, it is exactly the same. I also converted both the original and the new .kml file to a .svg using an online tool and compared them. They too were exactly the same.

Am I missing something here?


texts: Too long URL using visible layers list with arcgis js api WMSLayer

texts:
We have arcgis js api application (v 3.18) with TOC UI where users can turn individual layers on and off.

We have an issue with this WMS service. When users start turning on/off layers, WMSLayer sends a list of layers to WMS in URL parameter LAYERS=... The string with visible layers becomes too long and our app stops working on IE which has 2048 chars in URL limit, if I remember correctly. It still works on Chrome (it can handle longer URLs).

Is there a way to instruct WMSLayer to go switch to POST if URL is too long?
Do we need to have a proxy to handle POSTs as support for POST is not mandatory for WMS server? OR better is it possible to change URL to shorter one of GetMap GET requests and handle request on proxy communicating to real WMS using long URL?

Here is the example of failing REQUEST=GetMap URL. URL length is 2541. Chrome and Firefox open it normally: you see requested image. IE 11 truncates URL to 2048 chars, Edge truncates it to 2084 chars (strange number, maybe some Edge developer mistype the constant). So the shortened URL becomes invalid and you get error in IE and Edge.

Attempts to hook XHR communications with either

esriRequest.setRequestPreCallback(this.esriRequestPreCallback.bind(this));

or even brute force:

this.dojoXHR = dojo.xhr;
dojo.xhr = function (method, options) {
    console.log(method, options.url);
    return this.dojoXHR.apply(this.dojoXHR, arguments);
}.bind(this);

don't work. WMSLayer loads image needed by assigning URL to HTML Image element, i.e. image.src=url. These GETs don't go via XHR. I doubt if it's possible to hook to image.src setter.

Another approach to shorten URL is specification of GetMap params LAYERS. OpenGIS spec simply tells "The value of the LAYERS parameter is a comma-separated list of one or more valid layer names". ArcGIS WMS allows 2 forms: "Comma-separated list of one or more map layers (for example, 0,1,2 or layer1,layer2)." It's not clear how layer order is made, but this could help to communicate with this particular WMS which is running on ArcGIS. This GetMap request, for example uses short LAYERS form and it works.


texts: How to load Erdas .img files into Whitebox?

texts:
Just installed Whitebox 3.2.3 and I'd like to import a 14 MB.img file.

The data import options do not list .img as a file source.

A web search has also turned up nothing.


texts: Clear ArcGIS server rest cache with an in through script

texts: I would like to script the process to clear the rest cache on the \rest\admin page after I do a start/stop with the AGSSOM.exe provided on arcscripts.  Is it typical to have to clear the cache after creating new services or start/stop-ing services.  We seem to have to.  Thanks.


texts: Add Z and M values to stream network produced from DEM

texts:
I created a stream network from an existing DEM. My ultimate project goal is to be able to select a river from my stream network and use the 3D analyst to generate stream profile (Elevation by Distance).

The polylines I've created so far don't have M or Z values and I'm not sure how to add these. I have a DEM and I know I can add surface information which would add a z-max, z-min, and z-mean to the attribute table but that doesn't actually make it a polyline Z.

As of now to create a stream profile I only know how to use the 3D analyst by manually picking points/lines along the DEM. For my area this is very tedious and time consuming.

Any ideas?


texts: QGIS atlas-generation: Is it possible to preserve the map section?

texts: My aim is to generate several maps with the atlas function in QGIS 2.4 (win7 64) showing the same map section and scale. However, each generated atlas map shows a slightly different map section (not scale) according to the position of the features of the coverage layer. I want only the displayed features to change, not the map section. I can't see a possibility to 'fix' the map. Any ideas?


texts: What is a compatible data type for 'specific' date and time?

texts:
I need to import CSV file into PostgreSQL, but I got an error on time column which has data like this:

11/1/2016 12:00:00 AM

I tried to make column and assign 'timestamp' data type, but it's just not working.
Can anyone resolve this problem?


texts: How to copy information from Excel into QGIS Attribute-Table?

texts:
I have an Attribute-Table that I need to add a row from Excel into it.
I saved the Excel as CSV (UTF-8) and was able to copy and paste one slot at a time.

When I tried to copy and paste more than one slot nothing happens i.e.there was no error  message or anything like that.

I am sure that there is a much faster way apart from copying and pasting it one at a time.


texts: ArcGIS Pro - Rename attachments in attachment table based on relationship field

texts:
I'm looking for a Python 3/Aracade script to field calculate ATT_NAME in an attachment table based off of the REL_GLOBALID.
So where REL_GLOBALID is {67CECF26-9BE8-4E37-9408-3BC156DA51C0}, I'd like to rename ATT_NAME to Photo 1 and then where there is another attachment with the same REL_GLOBALID, I'd like that to be changed to Photo 2. See attached photo


texts: My Own Map server, Making my maps look like Mapbox streets or Skobbler

texts:
I want to make my own map server
However, I need to find the right style

The OSM style is pretty ugly, I really like the Skobbler maps http://maps.skobbler.com/ or Mapbox streets

I basically want to copy these designs and have my maps look the same
Are these styles publically available? does anyone have a link or method for obtaining them?

Is it as simple as downloading mapbox studio, selecting the Mapbox street template and exporting the styles and giving to the guy who is creating WMS ?


texts: Get all columns from postgis table when doing a getFeatureInfo request from mapserver's WMS layer

texts:
some background:

im using openlayers WMSGetFeatureInfo control to get feature info from a WMS layer.
the data is in postgis (its a table of points).
im using mapserver, and you need to define a template file for getfeatureinfo request.

if i put in the template file:

any string

it works. when clicking a feature i get back this string (for any feature of course).
if i put in the template file:

[date]

it returns the current date, (not a column in the table)
if i put in the template file:

[id]

it returns a number that is not the id column of the feature, though it is different for each feature.

how can i request all columns for that feature, as if i would do a "select *" quary ?
thanks allot!


texts: Error arcpy copy feature management problem when making join to shapefile

texts:
I am trying to join a csv file to a shapefile based on an ID field that keeps giving me this error:

ExecuteError: ERROR 000210: Cannot create output C:\Users\geo.gdb\2010reg
Failed to execute (CopyFeatures).

Here is my code:

import arcpy
env.workspace = r"C:\Users\geo.gdb"
arcpy.env.qualifiedFieldNames = False

# Set local variables for join
inFeatures = r'C:\Users\shapefile1.shp'
layerName = "ID"
joinTable = r"C:\Users\csv1.csv"
joinField = "ID"
outFeature = r"C:\Users\geo.gdb\successJoined.shp"

arcpy.MakeFeatureLayer_management (inFeatures,  layerName)

# Join the feature layer to a table
AddJoin_management(layerName, joinField, joinTable, joinField)

# Select desired features from veg_layer#
arcpy.SelectLayerByAttribute_management(layerName, "NEW_SELECTION")

# Copy the layer to a new permanent feature class
arcpy.CopyFeatures_management(layerName, outFeature)

I would like to copy all matched rows to a new output feature class. I do not want it to contain non-matching rows. The key "ID" is the same in both the csv and the shapefile. I am getting the error above that will not allow the final copy to work.


texts:
I have a database which i've exported as individual csv files for each table.

Below are some of the attributes which exist from two of the csv files:

wayPointID  linkID  order   x   y   z   length  radius

and

nodeID  x   y   z   name    description type

I'm in need of a method to convert this data into osm.

When I attempt just csv -> osm  in fme, of course the output is only the osm xml headers.

<?xml version="1.0" encoding="UTF-8"?>

I have been able to use 2DPointReplacer to take the x and y attributes and it successfully outputs nodes.

However, he data I have does require the tables be joined in order to build  attributes. Because of this I have to create a new id to each position and be ordered/continued across several tables.

I have written a python script which will do what I need but would like this done in FME.

Can anyone offer suggestions?

Thanks in advance.


texts: Viewshed Analysis for bridge polyline observer feature

texts:
I am using the Viewshed tool in ArcGIS 10.2 to determine the cells from a DEM that are visible from a highway corridor that has a bridge crossing a river (my observation feature).  My question is in regards to the z value of the observation feature, with several contingent follow-up questions:

does the viewshed tool interpolate the z value based on the DEM, or do I need to manually add a field and define this along the corridor?
If it does interpolate the z value, would I need to add OFFSET fields for the part of the corridor that is on a bridge structure (and therefore, not the same elevation as the DEM in that location)?
If I do need to add these fields, do I need to cut my line into multiple features/records to do this?


texts: What's the difference between Stereo 70 and Pulkovo stereo 70?

texts: I was wondering, can anyone tell me the diffrence between Pulkovo_1942_Adj_58_Stereo_70 (3844) and Dealul Piscului-Stereo 70 (31700)?


texts: How to extract features with specific coordinates from a ShapeFile?

texts:
I would like to extract a number of features from a ShapeFile which match with specific latitude/longitude values.  I am interested if/how this can be done with QGIS and ogr2ogr. The output should be a ShapeFile as well; for QGIS a layer is sufficient.

I assume the command for ogr2ogr could be something like:

ogr2ogr -f "ESRI Shapefile" -where "Y = 40.752459 AND X = -73.983521" output.shp input.shp

However I am having difficulties finding out what the attributes for lat/lon are named.

Optional: Also, is there an option for a fuzzy query where I can configure a radius around the coordinates?

I am running Ubuntu 12.04. with QGIS 1.9.0 via ubuntugis-nightly and GDAL 1.10.0.


texts: Problem updating Georeferencing for TIFF

texts:
When I update georeferencing for TIFFs, the map disappears and doesn't save. Auto adjust is off, rectify wont save, control points I can save seperately. Settings are set to documentation suggestions. How can I update and save these images? My workflow goes as follows:

Add data(tiff file), asked if I want to create pyramids to which I agree, fit to display, use 2-3 control points to line up, then I usually update processing, but I tried to rectify and it pops up saying failed to save even though the .ovr file has been created and in the same location as original TIFF. I was on the phone with Esri yesterday who said all I needed to do was turn off auto-fit (which I did and update geoprocessing actually worked!) but now today I am back at the same issue. I have also restarted computer/cleared cache several times


texts: How to export only labels as svg in QGIS 2.6?

texts:
I'm working with a big big data set in QGIS 2.6 (but also in QGIS 2.7).

I have many line layers styled with the rule rendering tool. What I want to do is to export only the labels of all these line layers as svg while the line style can be exported as png.

I tried a small work-around: I have set a total transparency for the line layer leaving only the labels. The problem is that during the svg exportation, QGIS seems to read all the line anyway and these causes an error (I suppose because the dimension of the whole file).

Is there a way to do what I need?


texts: See python script running while processing a task in QGIS?

texts: Is it possible to see scripts running in Python Console - everything from tools and analysis vector to raster?


texts: InSAR processing: why can't we just compute distances directly from time delays?

texts:
My question can be very basic, but even after studying quite a few tutorials on InSAR processing, I still do not understand one thing: if we know a wavelength, speed of light, a time delay it takes the signal to travel to each ground cell and back, and even precise locations of the satellite - why can't we just compute distances for each pixel?

In other words, why do all those interferograms are involved in computation of displacements and topography, if all we need is (given all the mentioned variables) just to calculate distance(s) and some trigonometry?


texts:
R is becoming pretty strong tool for handling and analyzing spatial data. I learned some useful things through questions like these at SO and thought it might be useful to have something simmilar, but more 'spatially' oriented.

Can you share some spatial R tips and tricks that you found useful?


texts:
How accurate is the Open-Elevation API?
I'm trying to follow this example (except with s=1000 to get high-accuracy data), and get the following elevation profile:

However, using Google Earth Pro, I get a much finer elevation profile (the accuracy of the start and end points appear OK):

Do we know why that is? Does Google use a different dataset, or does it essentially apply some interpolation to the results I get with the open-elevation API? The location is Greece in case that matters.
I'm also unable to set s=2000 or higher for some reason, so I can't increase the resolution more apparently:
Traceback (most recent call last):
  File "/Users/username/Desktop/Radio Astronomy/linkCoverage/example_elevation.py", line 64, in <module>
    fp=urllib.request.urlopen(response)
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/urllib/request.py", line 222, in urlopen
    return opener.open(url, data, timeout)
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/urllib/request.py", line 531, in open
    response = meth(req, response)
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/urllib/request.py", line 641, in http_response
    'http', request, response, code, msg, hdrs)
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/urllib/request.py", line 569, in error
    return self._call_chain(*args)
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/urllib/request.py", line 503, in _call_chain
    result = func(*args)
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/urllib/request.py", line 649, in http_error_default
    raise HTTPError(req.full_url, code, msg, hdrs, fp)
urllib.error.HTTPError: HTTP Error 400: Bad Request


texts:
In QGIS 3.0.x, there was a menu

Project -> Save Map As Image

which allows one to export a portion of the current view (screen shot) into a png file at a given dpi.

It is quite convenient, but it seems to be missing in QGIS 3.2.0-Bonn (revision bc43194)

Is it moved to somewhere else from the Project menu or removed/disabled/hidden?

How do I get the equivalent function or re-enable it?


texts: Does anyone know of a way that you can style OS Zoomestack and order the layers so that each time it is inserted in QGIS 3 it styles automatically and in the correct order?


texts: Unable to pan on map on a specific zoom level . Openlayers2

texts: I am fixing bug of a big application, In which after I cant pan over map but after zoom in I can pan. after again zoom out Its impossible. As I mentioned its big application so I cant post all code here . Any one have idea where error can be ? Before negative vote kindly tell me how I can improve question as I dont know where error could be ?


texts:
Is it possible to use "and" & "or" queries on a WFS service through the URL?

Here's my current URL:

http://mywebsite.com/county_roads_wms_only/MapServer/WFSServer?request=GetFeature&service=WFS&typename=county_roads_wms_only:Roads&maxfeatures=2&filter=<ogc:Filter><ogc:PropertyIsLike escape="\" singleChar="_" wildCard="%"><ogc:PropertyName>Road_Label</ogc:PropertyName><ogc:Literal>%mcculloch%</ogc:Literal></ogc:PropertyIsLike></ogc:Filter>'

Is there a way to form it so that it would be "%mcculloch% AND %st%" or "%mcculloch% OR %st% in the URL?


texts: QGIS, Tilemap stored in local folder or in project directory

texts:
I created an offline tilemap and want to use it as background map in QGIS. I like to store the maptiles in the project folder. The QGIS-Project should work independently from any installation or plugin.

How can I get this to work?


texts: I want to create NDVI of multiple raster using ArcPy

texts:
I have multiple raster data which are layer stack. Now I want to create NDVI of each raster data. I am stuck in the code, I developed a script but dose not generate any output.
import arcpy, os
from arcpy.sa import *

arcpy.env.workspace = "L:\\Arcpy\Data\\raster_data"
arcpy.env.overwriteOutput = True

input_raster = arcpy.ListRasters()
outputDir = ("L:\\Arcpy\\Data\\NDVI")
print(input_raster)

for data in input_raster:
    red = arcpy.MakeRasterLayer_management(data,"L:\\Arcpy\\Data\\raster_data\\Anugul_red.img", "", "", "3")  # Name changed
    nir = arcpy.MakeRasterLayer_management(data,"L:\\Arcpy\\Data\\raster_data\\Anugul_nir.img", "", "", "4")  # Name changed
    # the first part of the file name before extension (if present) then add .img
    output = os.path.join(outputDir,os.path.splitext(data)[0].split('_')[0] + '_ndvi' + '.img')
    ndvi = (Float(nir) - Float(red)) / (Float(nir) + Float(red))
    print("{}\n{}".format(output, ndvi))
    ndvi.save(output)


texts: How to replace the leaflet-draw toolbar by buttons outside the mapDiv

texts:
I did quit research to move the bar out of the leaflet map.
So what I did basically is that I created 3 buttons outside the map Div:

Draw Polygon
Delete
Edit

I binded events to the draw Polygon with :

new L.Draw.Polygon(this.map, this.drawOptions.draw.polygon).enable();

But I couldn't do the same with Delete and Edit. I found in a forum that I could use:

new L.EditToolbar.Edit(map, {
            featureGroup: drawControl.options.featureGroup,
            selectedPathOptions: drawControl.options.edit.selectedPathOptions
        })

But I didn't find the typings for L.EditToolbar since I'm developing in Ts. Any Ideas ? I want also to remove the Toolbar from MapDiv without touching to the draw options since they are used by my Draw buttons.


texts:
I have local Overpass server installed with some snapshot data.

Is there something like query data into graphical editor, edit the geometries and then update the Overpass server with the modified data?


texts: Converting Esri Shapefile to Geojson file

texts: I am trying to convert an Esri shapefile, e.g. wards boundary data, to a Geojson file format.


texts: Ground Penetrating Radar (GPR) to ESRI raster format

texts:
can anyone recommend a software package to convert ground penetrating radar data (in .dzt format) into an ESRI raster format? - Preferably free software.

I have come across ZondST2D, but as far as I can tell it will not do the conversion.


texts: Copying coordinate system and extent from a TIFF to a polygon shapefile in SAGA GIS

texts:
I used R to batch process object based image segmentation on some geoTIFF files using rgdal to communicate with SAGA GIS. I got the resulting shapefiles, all with the following properties:

Projection Undefined Coordinate System
West -0.5
East 7999.5
West-East 8000
South -0.5
North 3999.5
South-North 4000
Type Polygon
Vertex Type X, Y

For each individual file, I got the following warnings:

15:13:13: Warning: Unknown field with tag 33550 (0x830e) encountered (in module "TIFFReadDirectory")
15:13:13: Warning: Unknown field with tag 33922 (0x8482) encountered (in module "TIFFReadDirectory")
15:13:13: Warning: Unknown field with tag 34735 (0x87af) encountered (in module "TIFFReadDirectory")
15:13:13: Warning: Unknown field with tag 33550 (0x830e) encountered (in module "TIFFReadDirectory")
15:13:13: Warning: Unknown field with tag 33922 (0x8482) encountered (in module "TIFFReadDirectory")
15:13:13: Warning: Unknown field with tag 34735 (0x87af) encountered (in module "TIFFReadDirectory")

The shapefiles that were created are fine, expect that when I attempt to project them in any given coordinate system, it projects incorrectly. I am not sure what the standard or default projection is in SAGA, as it appears that when the TIFFs were imported, they lost their spatial reference/coordinate system.  I have tried defining and projecting a coordinate system, but when I change the extent of the grid file, I get an empty map.

I would like to copy the coordinate system and extent of the TIFF file and "assign" or coerce it to its respective shapefile in SAGA. The desired attributes are as follows:

Projection
Projected Coordinate System [EPSG 2006]:
St.Lucia1955/BritishWestIndiesGrid
[+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +towgs84=-149,128,296,0,0,0,0 +units=m +no_defs ]
West
499993.1
East
502004.5
West-East
2011.400000000023
South
1522994.7
North
1524000.4
South-North
1005.6999999999534
Cell Size
8.9


texts: Extracting area from map and making this new layer?

texts:
It is probably the most simplest thing to most people, but I am only starting out in the GIS world. I have a map of Ireland but my focus is on one particular county. I want to 'crop' the rest of Ireland out and just view my county of choice.

How can I go about doing this?


texts: How to create a reservoir for a dam site?

texts:
I have created the watershed for a proposed dam site point using SRTM 90m DEM. Now my task is to create a reservoir or pondage where water would be collected from the watershed.

Above a watershed can be seen with a red line on which Dam can be created. A dam also has a reservoir where water is collected and it is always inside the watershed, now problem is to delineate that reservoir.

Practically, reservoir should be in lowest elevation (flat area), so can we create such using Dem data.


texts: Strategies and methods that Datums and Coordinate Systems use to account for movement of the earth's surface

texts:
Are there datums and coordinate systems that address the issue of shifting coordinates through time? If so, what are some of the strategies and methods that the datums and coordinate systems use to address the issue?

To provide some context:

Several Earth processes, plate tectonics and isostatic rebound, affect where a (x, y, z) location is (or perhaps better, is calculated to be) in a time step t1 in comparison to a different time step t2. In some areas of the world, such as California, USA, this movement can be in the magnitude of centimeters per year. For some GIS analysis the shift in (x, y, z) locations could have a serious effect on the results, for example examining sediment flux that rely on comparing raster cells and/or coordinates across different time steps often spanning several years or decades.


texts: How to run in background an .exe program withArcGIS?

texts:
I am working with the ArcGIS program that reads Python language. For that I've already implement scripts as script tool on the ArcGIS.
Now I need to run in background an exe program "gams".
Is it possible to do directly from the ArcGIS? Or is it possible to do with a new python script? And how can I do this?


texts: Listing multiple fields to be dropped at once in QGIS modeler?

texts:
I'm trying to drop more than one fields at once in QGIS 3 modeler on Microsoft Windows. I'm using the algorithm "Drop field(s)".

What is the good syntax to list multiple fields to be dropped?

Not working:
'Field1','Field2','Field3'

"Field1","Field2","Field3"

When I run it, I get a python error and the result is empty, no column at all.


texts: Loading a raster using QGIS into a PostGIS2.0 enabled database

texts:
Following on from my attempt to install a postgresql 9.1 db with postgis 2.0 on windows 7, I am trying to figure out how to load a raster. I have successfully loaded a shapefile, and am trying to do the same with a .tif file using the Load Raster to PostGIS plugin (version 0.5.1) in QGIS (version 1.7). I have set up the connection to my db, and am using the following settings:

When I click on the OK button, I get the below error message. I've tried this with a .adf file and also a .tif file, both projected in Albert Equal Area Conic with an SRID of 102003.

Checking parameters...
Connecting to database...
Storing overview 1 on database...
Failed.
Finished storing overview 1.
Storing overview 2 on database...
Failed.
Finished storing overview 2.
Storing overview 3 on database...
Failed.
Finished storing overview 3.
Finished.

This process has not inserted anything into my database, and I don't understand the error message. Some previous related questions are here (asked by me, using a different process on loading rasters) and here (using the same plugin but much earlier in the year).


texts: How to find Shortest Route between two points in PostGIS without using pgRouting?

texts:
I want to find shortest route between two points stored in PostGIS database. I know about pgRouting functionalities but I don't want to use it because of data processing/topology creation processes.

Is their any such function/functions in PostGIS to achieve the same?


texts: What would cause a Cesium 3D-Tiles tileset to not even request it's content?

texts:
I am generating an I3DM through Python. I know that the I3DM is probably ok, since the 3d-tiles-tools allow me to extract the binary GLTF file from it. However, when I go to render the tileset in Cesium, it shows in the Web Tools Inspector it doesn't even request the i3dm file. It requests the tileset.json, but that's all. What would cause this to happen? Here's my tileset:

{
  "asset": {
    "version": "1.0"
  },
  "properties": {
  },
  "geometricError": 70,
  "root": {
    "refine": "ADD",
    "boundingVolume": {
      "region": [
        -2591020.1425219355,
        1588193.8474433378,
        318979.8574780645,
        -3805806.1525566624,
        0,
        100000
      ]
    },
    "geometricError": 0,
    "content": {
      "uri": "myi3DM.i3dm"
    }
  }
}

The I3DM's feature table only has the POSITION variable in it's binary.


texts: Prevent icons from flashing during ArcGIS FeatureService layer refresh

texts:
Here is my current situation:

Web Mapping app built with ESRI WebAppBuilder 2.10
FeatureService of geometry type Point with Feature Access and Web Mapping capabilities enabled generates features from table on PostGIS database. Uses the column PK as the ObjectID. Contains ~ 400 features on average that may have any of their properties changing except for ID.
Web Map contains FeatureService and draws features markers on map.
During layer refresh on 5 second interval, all the features are redrawn, the icons blink off and then back on, giving the app a very janky feel.

Is it possible to prevent the features from removed and then re-drawn to the map during layer refresh? I couldn't find anything specific to this in ESRI documentation, and the map layer drawing logic is contained by default within ESRI's CDN, and is difficult to access through the file available in the WebAppBuilder application.

I've tried disabling custom symbology (I've assigned custom icons within the layer definition) and using different unique values for the ObjectID, but still no dice. It seems like maintaining a live FeatureService that is updated frequently is a fairly basic use case, but I can't find any kind of support for this on the ESRI forums.

Any ideas?


texts: PostGIS ST_Intersects does not return expected result

texts:
I want to search objects that have at least one point inside of a rectangle. Based on the size of the rectangle, PostGIS sometimes returns the contained object, sometimes not.
create table t1 (id integer not null primary key, coords geography(linestring));
insert into t1 values (1, 'linestring(9.4 48.5, 9.2 48.6, 9.0 48.8)')

So the table contains a linestring, that consists of 3 coordinates (somewhere in the area around Stuttgart).
Now I want to search all objects, that have at least one point inside of a given rectangle. I use a polygon for this boundary. The following query returns the linestring from above:
select t.id, st_astext(t.coords) from t1 t
where ST_Intersects(
    ST_GeogFromText('POLYGON((-110.0 49.0, 10.0 49.0, 10.0 47.0, -110.0 47.0, -110.0 49.0))'),
    t.coords)

> 1, LINESTRING(9.4 48.5,9.2 48.6,9 48.8)

In the following query, I just modified the right (eastern) border of the rectangle and moved it from lon=10.0 to lon=20.0. This query does not return any result:
select t.id, st_astext(t.coords) from t1 t
where ST_Intersects(
    ST_GeogFromText('POLYGON((-110.0 49.0, 20.0 49.0, 20.0 47.0, -110.0 47.0, -110.0 49.0))'),
    t.coords)

What am I doing wrong?
Some info about my environment:
SELECT version() || ' ' || postgis_full_version();
PostgreSQL 12.3 (Debian 12.3-1.pgdg100+1) on x86_64-pc-linux-gnu, compiled by gcc (Debian 8.3.0-6) 8.3.0, 64-bit POSTGIS="3.0.1 ec2a9aa" [EXTENSION] PGSQL="120" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 5.2.0, September 15th, 2018" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.3.1" WAGYU="0.4.3 (Internal)"


texts: HIFLD Railroads DB contains coordinates that don't make sense

texts:
I downloaded the HIFLD Railroads DB and imported it into MySQL with OSGeo4W. I then do this query and the coordinates I get just don't make sense:
SELECT ST_AsGeoJson(SHAPE)
FROM railroads
WHERE OGR_FID = 123751;

Here's what I get back:
{
   "type":"LineString",
   "coordinates":[
      [
         -10791447.5365,
         3527298.6227
      ],
...
   ]
}

That's not a GPS coordinate. idk what that is.
When I use OSGeo4W with another DB (TIGER County DB) I get coordinates that do make sense. The query I use with that:
SELECT ST_AsGeoJson(SHAPE)
FROM county_shapes
WHERE OGR_FID = 1637;

The result:
{
   "type":"Polygon",
   "coordinates":[
      [
         [
            -97.997996,
            30.720209
         ],
...
      ]
   ]
}

Those coordinates make sense so using OSGeo4W can work but it just isn't for the HIFLD Railroads DB.
Any ideas?


texts: Read date of file creation from within a zip file using FME

texts:
I have a simple process that downloads a file from an FTP site and unzips the contents to a folder. Unfortunately the date of the file once unzipped is recorded as the date when the fmw was run but I need to get the date of the file itself. Is this possible?
The process

The result

But looking at the zip file we see the date that we need.

I noticed the zip extractor has a _creation_instance can this be captured into a variable (assuming it is the date we need)?
Is there a parameter in the ZipExtractor (or other plugins) that allows me to capture the date of the files within it?
UPDATE:
Based on the responses I tried the following but the powershell doesn't run from FME and without it I get the wrong date values.


texts: I took an aerial image of our farm using RC quadcopter. How am I going to orthorectify this image? Please can you help me with detailed instructions/methods involving the process.


texts: Features don't appear importing shapefile into QGIS

texts:
I am trying to import a shapefile into my map, but nothing appears on the map. When  I look at the attribute table the top line reads: "filename:: Features Total: 117412, Filtered: 0, Selected:0" and all of the expected fields are there, but there are no features visible in the attribute table either.

In the attribute table I have selected Show All Features, on the map I have selected Settings/Options/Data Sources -> Attribute table behaviour = Show all features. When I right click on the layer in the Layers tab and say Zoom to Layer the map zooms to approximately what I would expect the extents of the layer to be.

Overall it looks like the features are all there, but somehow they are not visible in either the attribute table or on the map. I actually received two shapefiles from the same source at the same time, the first is displaying as I expected it to, it is only the second of them that is giving this problem.

I'm running 3.8.3 - Zanzibar (although I thought I had updated to the latest version just the other week. I'll do that now though, seeing that 3.12.0 is available).
Any suggestions would be greatly appreciated.


texts: How do I make sure when importing x and y data that they are in the same coordinate system as my basemap?

texts:
I have an ArcMap basemap of Britain which has a coordinate system which is National Grid. I am trying to import a spreadsheet of coordinates of sites across Britain which are in decimal degrees (eg 50.933392 -1.5033064). However when I use the convert to x and y tool my data points are all projected into the sea off Cornwall. (Presumably my coordinate system is wrong so it puts the point in metres from a 0, 0 point.
How do I sort this out?


texts: SearchCursor stops after first iteration

texts:
I'm working on a script that takes a few feature classes, unions them together and then runs a series of geoprocesses on them and puts it into a single feature class. From that feature class I have a for loop set up that is to use CopyFeatures_management to create new feature classes based on the value of a specific field (in this case: PIN) and name the new feature class "CAUV_"+PIN(from attribute table).
It runs through everything once and creates the new feature class as it should but then it stops. It won't iterate through the rest of them. What on earth did I miss?
#Iterate over values in attribute table
def unique_values(table, field):
    cursor = arcpy.da.SearchCursor(table, [field])
    unique_values = []
    for row in cursor:
        if not row[0] in unique_values:
            unique_values.append(row[0])
            return sorted(unique_values)
unique_pin = unique_values(table=ParcelsLandUse, field="ParcelPublishingAuditor_PIN")

#Select Features and Make Layer
ParcelsLandUse_lyr = arcpy.MakeFeatureLayer_management(ParcelsLandUse, "ParcelsLandUse_lyr")
with arcpy.da.SearchCursor(ParcelsLandUse_lyr, "*") as cursor:
    for PIN in unique_pin:
        expression = "%s = '%s'" % ("ParcelPublishingAuditor_PIN", PIN)
        arcpy.SelectLayerByAttribute_management(ParcelsLandUse_lyr, "NEW_SELECTION", expression)
        outCopy = "Y:\\Data\\AuditorGIS.gdb\\Calculations\\CAUV_"+PIN
        arcpy.CopyFeatures_management(ParcelsLandUse_lyr, outCopy)
        arcpy.SelectLayerByAttribute_management(ParcelsLandUse_lyr, "CLEAR_SELECTION")
        info = "Processed: " + outCopy
del unique_pin
del PIN
info = "CAUV Parcel Soils Complete"
arcpy.AddMessage(info)


texts:
I made a custom map tool plugin with pygis like this:

class myTool:
    def __init__(self, iface):
        self.iface = iface
        self.canvas = self.iface.mapCanvas()
        ...

    def run(self):
        self.mapTool = mapTool(self.canvas)
        ....

class mapTool(QgsMapToolEmitPoint)
    def __init__(self, canvas):
        self.canvas = canvas
        QgsMapToolEmitPoint.__init__(self, self.canvas)
        self.rubberband = QgsRubberBand(self.canvas, QGis.Polygon )
        ....

And i make some calculations inside the mapTool class and i want to display the results into Qgis status bar, but i haven't found a way to accsess the status bar from the mapTool class because it uses the canvas attribute not the iface attribute.


texts: Encode MapServer tile pixels with data values

texts:
Is it possible to encode values from raster bands into pixel output using a MapServer .map file?

More and more I'm finding that a commonly adopted method for moving geographical data for rendering WebGL based maps is the old game developers trick of using PNG/JPG images with float/integer values encoded into the RGBA bytes. These images are then fed into GLSL shaders as textures where they are decoded to the original (but slightly offset) values on the GPU.

There are various projects using this method,and subsequently various utilities that are being used to encode these values. e.g. https://blog.mapbox.com/how-i-built-a-wind-map-with-webgl-b63022b5537f

An example of an encoding algorithm is:

png.r = Math.floor(255 * (u[i] - u.minimum) / (u.maximum - u.minimum));
png.g = Math.floor(255 * (v[i] - v.minimum) / (v.maximum - v.minimum));

This would be an incredibly convenient way of passing texture data directly to GLSL shaders using common mapping tools such as Leaflet, OpenLayers or Mapbox, without the need to pre-process image files

I can generate a greyscale image using

  LAYER
    NAME "encode" STATUS ON
    DATA "vstm-uv.grb2"
    PROCESSING "BANDS=1,2"
    TYPE RASTER
    STATUS ON
    PROJECTION
      "+proj=longlat +a=6371229 +b=6371229 +no_defs +lon_wrap=180"
    END
  END #Layer

But i can't figure out how to color these values or process the incoming [red] & [green] bands.


texts: Adding attributes of points to attribut table of lines in QGIS

texts:
I'm using QGIS 3.8 with GRASS 7.6.1
I want to add only the names of my vertices (from line intersection, duplicates deleted) to the attributtable of my lines (like start and end node). Every line has only two vertices.
What I tried:

geom_to_wkt(start_point($geometry)) --> gives me the coordinates (I need the name) of my real nodes not of the vertices

"add by position" only transfers the name of one vertices (but I need both)


texts: How to write FIDs of points within each polygon to polygon field?

texts:
I have two feature classes one with points one with polygons (in ArcGis 10.2). Now, I want to write into a String Field of the polygon layer a list of all the Point-FIDs that are within  each polygon (often more than one, which is why join doesn't help).
I was hoping that there is a possibility using Field Calculator and somehow concatenate the FIDs. But I don't know how to access the other feature class within the calculation code - is that possible? Also a Python code would be appreciated or anything else...


texts: Change the starting measure point of a route - linear referencing

texts: I have a few routes created and stored as a feature class in a GDB.  When using the create route tool there is an option to choose the starting point of the routes.  However, I only have certain routes that start TOP LEFT (of the bounding rectangle), some start LOWER LEFT etc etc.  Is there a way to change the starting measure point after you have created routes or am I going to have to create routes in groups (those with lower left start points, the upperleft etc) and then merge afterwards?


texts: Telling Tabulate Area to use custom textual labels for column names?

texts:
I have a land use raster with unique values of 1-8 and the appropriate label names for the land use are defined under the symbology tab.

When I process that raster with the Tabulate Area tool in conjunction with census divisions, the column names of the resulting table are VALUE_1, VALUE_2, ..., VALUE_8. Is there a way to create such a table but with the proper textual labels for each type of land use instead?


texts: Query out a specific date regardless of the time, using python

texts:
Problem: I have a feature class with an attribute table with a field called Edited_Date and its data in the field is as follows: 4/8/2015 10:00:00 PM. I want to select all records that are Edited_Date is 4/8/2015 regardless of the time of day. SO I want to use a wildcard to substitute the hour minutes seconds and PM.

I run the following code using python to select records by Edited_Date but no records are being selected. I've tried different formats but none is working and I can't find the proper syntax to query pout a date such as

Code:

arcpy.SelectLayerByAttribute_management("Plant_Layer", "NEW_SELECTION", " [last_edited_date] = '4/8/2015' + '%s' ")

What is the syntax to do this in Python? Select a date regardless of the time.


texts: Adding up several rasters to show the change

texts:
I want to add several images such that their change is shown in different colors. Such as the following in which I simply added two  images in their grayscale. But I would like to have a colored representation of change when several images are added together. For instance for 50 images, I have their difference depicted in 50 different colors.

My code to add grayscale images is as follows:

imageA = cv2.imread('/2013.tif')
imageB = cv2.imread('/2014.tif')

grayA = cv2.cvtColor(imageA, cv2.COLOR_BGR2GRAY)
grayB = cv2.cvtColor(imageB, cv2.COLOR_BGR2GRAY)

result = cv2.add(imageA,imageB)

Is there any opensource approach using Python/R for that?


texts:
I am trying to build a web app where users can use OpenLayers to draw features and then save them either as a download or write to a file on the server which I can collect when I need it. I am happy with the OL3 bit so far.
From reading around it seems AJAX is the way to do this so I have managed to included a button in my HMTL which runs a function that runs a php file to create a new file.I was pretty happy that I managed to get an empty text file as a start. So how do I get that file to contain the JSON information? I assume I need to send it from the JS world to PHP somehow?

The Function

function loadXMLDoc()
{
var xmlhttp=new XMLHttpRequest();

xmlhttp.open("GET","run.php",true);
xmlhttp.send();
}
</script>

The button

<button type="button" onclick="loadXMLDoc()">do something</button>

The PHP file (run.php)

<?php
$myfile = fopen("store/anothernewfile.geojson", "w") or die("Unable to open file!");
$txt = ['lyr_site'];
fwrite($myfile, $txt);
fclose($myfile);
?>

I don't have a geoserver and don't think it is really necessary for something simple like this, using .geojson files seems to be working so far.


texts: What is the history behind WKT and proj.4?

texts: What is the history behind WKT and proj.4 formats for defining coordinate systems? Does one have advantages over the other? Why do both exist? I'm guessing they were developed independently for separate systems, but I'd be curious to know more specific reasons.


texts: Example of C++ code in Visual Studio initializing a created raster

texts:
I refer to Sample for MS Visual Studio C++ using GDAL library

I noticed that the answer mentions that when GDAL is built from source there would be a tools folder with examples of how to initialize a created raster. I have GDAL built from source (following the steps specified here: http://crowdsimulation.blogspot.com/2015/07/geospatial-data-abstraction-library.html?showComment=1452035822922#c1369732198978557869) in Visual Studio 2010 and a 32 bit machine.

However, I am not able to locate the tools folder with the examples.
Is there any other source where I can obtain the examples?
Can someone post a simple one?


texts: Compute total area of several polygons within a polygon

texts:
I have joined a dataset containing the polygons of green spaces, with a dataset containing vegetation polygons using intersect. The result is a dataset with separate entries for each vegetation polygon within the green spaces. Looking like this:

What I want is the total area of each vegetation type within each of the green spaces, similar to this (this table came from a trial using the tabulate areas, right format but wrong areas):

What would be the best way to achieve that, using python? The area would have to be calculated for each vegetation type (GRIDCODE) within each Green space. The areas are then meant to be attached to my original dataset containing a lot more information about the green spaces.
I've tried using a loop to compute it separately for each row, but did not get it to work.


texts: Wrong coordinates but polygons and map are aligned just fine in QGIS

texts:
I have 2 polygons I downloaded from public databases and I used the Bing aerial map with the OpenLayers plugin. They are aligned PERFECTLY, they are in the area they are supposed to and everything was great. Everything is in EPSG:4326 - WGS 84. The problem is when I did the composition I needed the coordinates to show up and... they are wrong, they simply do not reflect the real coordinates.
For reference, I took from literature the coordinates of the centroid of my principal polygon:  25°29'28.4350", -100°31'38.9740", also I put this coordinates on Google Maps and they are correct. But my map (composition) show completely different coordinates 6°51'0.000"  11°0'0.00" (no exactly in the centroid, but they still very different from each other)...
I think I'm using the same format of coordinates (WGS 84) but I am very confused, I don't know how to change the coordinate in order for they to match.
Also, I obtained the coordinates from the centroid of my polygon (to see what is the problem) in the properties table and obtained coordinates that are way closer but still off: y=25 x=-101.
I just have no clue to what to do.


texts: Converting Between Similar Classes in ArcObjects and WPF Runtime Libraries (C#)?

texts:
There seems to be an "impedance mismatch" (if I can borrow the term) between similar classes in ArcObjects and the WPF runtime libraries. For example, there is a ESRI.ArcGIS.Client.Geometry.Polyline in the WPF runtime and a ESRI.ArcGIS.Geometry.IPolyline (and corresponding PolylineClass) in ArcObjects.

I want to be able to do something like this (or at least similar):

var poly = new ESRI.ArcGIS.Client.Geometry.Polyline();
var ipoly = new ESRI.ArcGIS.Geometry.PolylineClass();

// stuff happens

poly = ipoly;

Is there an easy way to convert between the two? Or am I going to have to write some deep copy methods?


texts: QGIS modeler: using input rasters filename for naming columns by expression

texts:
I am running QGIS 3.12.3 București on Windows 10 x64; I am building a simple model to reclassify a series of rasters with a 2x3 table, and then count each unique reclassified value within different polygons in a shapefile.

My input rasters follow a naming convention with area_type_year_satellite, and I have rasters for each year over a long time span.

Now, since the input raster for the zonal histogram has to be a single entry, I am running the model as a batch process. I am successfully saving the output shapefiles with unique names according to the year in the filename of each input rasters, as suggested here, using an expression similar to:

'C:/Users/user/Desktop/' ||  @inputsitespolygon || '_' || substr(@RastertoReclassify, 10, 4) || '.shp'

My issue is: the Zonal Histogram algorithm allows for predefining the name of each column of the output shapefile, and I am trying to obtain unique column names depending on the input raster filename, with a similar approach to the above, using either:

substr(@RastertoReclassify, 10, 4) || '_'

or

concat(substr(@RastertoReclassify, 10, 4), '_')

However, any expression I used seems to fail, and the columns in the output keep the default name (HISTO_) or register only the last part of the concat function:

My question is: is there a way to programmatically name the output columns based on the input raster filename?

Ideally, they should be something like  1986_0 | 1986_1 for one output; 1987_0 | 1987_1 for another, and so on. I have to say I am not familiar with Python, and I don't know how to build a custom processing script. Not sure if I am just missing something, but I am a bit clueless on how to achieve this.


texts:
I'm using QGIS 3.4 , and without going into python scripting, is there a way I can create some kind of rule based styles where I can add a parameter?

I have a number of layers, each with between 1 and 10 rule based styles. For instance, here is one with 3 rules:

Plan1818 = 'Y' AND ("Type" = 'Deposit') and Status = 'Used'
Plan1818 = 'Y' AND ("Type" = 'Cut') and Status = 'Used'
Plan1818 = 'Y' AND ("Type" = 'Structure') and Status = 'Used'

I have 20+ 'Plans', and growing, and I am loading the saved style from one, loading it the next, and then changing the plan number for each rule.

I feel as there should be way to add the plan number as a parameter, but I've not found one yet. I tried a separate, rule Plan1818 = 'Y' and removing it from each of the others, but that didn't seem to work, even adjusting the symbol level.

Suggestions?


texts: Missing multipolygon after ImpOSM import to PostGIS

texts:
I've been testing ImpOSM 2.5.0 with an extract of Washington State (in the U.S.) and am very impressed with the overall result.  However, I noticed that fairly prominent lake was missing from the resulting "osm" DB, so I started to dig a bit deeper.  I confirmed that the lake is present in the extract, here are it's details:

    <relation id="411923" version="13">
            <member type="way" ref="42114354" role="inner"/>
            <member type="way" ref="42114348" role="inner"/>
            <member type="way" ref="42114344" role="inner"/>
            <member type="way" ref="42114346" role="inner"/>
            <member type="way" ref="42115602" role="inner"/>
            <member type="way" ref="42117265" role="inner"/>
            <member type="way" ref="42125939" role="inner"/>
            <member type="way" ref="52374472" role="outer"/>
            <member type="way" ref="126149886" role="outer"/>
            <member type="way" ref="126149894" role="outer"/>
            <member type="way" ref="52434337" role="outer"/>
            <member type="way" ref="126263156" role="inner"/>
            <tag k="name" v="Moses Lake"/>
            <tag k="natural" v="water"/>
            <tag k="source" v="yahoo_wms"/>
            <tag k="type" v="multipolygon"/>
    </relation>

I then checked the DB table definition in geometryc_columns and saw this:

  osm=# select * from geometry_columns where f_table_name like '%water%';
   f_table_catalog | f_table_schema |      f_table_name       | f_geometry_column | coord_dimension |  srid  |    type
  -----------------+----------------+-------------------------+-------------------+-----------------+--------+------------
                   | public         | osm_new_waterways       | geometry          |               2 | 900913 | LINESTRING
                   | public         | osm_new_waterareas      | geometry          |               2 | 900913 | GEOMETRY
                   | public         | osm_new_waterareas_gen1 | geometry          |               2 | 900913 | GEOMETRY
                   | public         | osm_new_waterareas_gen0 | geometry          |               2 | 900913 | GEOMETRY

The GEOMETRY type is new to me, but I assume that it is a valid type to store polygons/multi-polygons.

I followed the tutorial at http://imposm.org/docs/imposm/2.5.0/tutorial.html and had no problems with any of the steps listed there, but I'm curious if given the nature of this feature, should I have done something more to enable it to be parsed from the extract and created in the DB?

Many thanks for any clues you can provide.


texts: Open Source Desktop Tool to Style MapBox Vector Tile Data for MapBox

texts:
I have MapBox Vector Tile data in the form of mbtile files that needs to be styled for MapBox.  Need a desktop application to work with data locally that will allow the style to be saved locally as style.json.  I have looked at QGIS, but it does not appear to let you save as styles compatible with MapBox.
Any ideas on how to accomplish this task?


texts: How do I convert a shapefile feature class from Angular Unit: Degree to linear feet?

texts:
I have some ground receptors that are in decimal degrees and I need them to show linear feet instead. Does anyone know how to go about this?

GCS: GCS_WGS_1984
Datum: D_WGS_1984

Also, can someone please create a tag for LinearFeet?


texts: Creating multiple random points per polygon based on attribute table in QGIS?

texts:
I have read Creating random points in QGIS with values from Attribute Table, but it doesn't address my question because I use the same process:

Vector > Research Tools > Random points inside polygons (variable)

I am presented with the following console:

The attribute NUMFARMS is in the attribute table, is integer-valued, and contains 0s for most of the features. Some of the features have a value higher than 1.

The randomly generated points, however, are only ever one per polygon. What am I doing wrong?

EDIT: In response to a comment, here is a screenshot of the Fields tab:


texts: Why do I get this anomaly (vertical line) in my TWI raster?

texts: I have derived TWI using SAGA GIS, Have a look on the center of the image, a vertical line goes through it from top to bottom. Anyone else confronted this issue?


texts:
I have a QGIS 3.2 workflow that helps me to mass produce image classiication labels for some deep learning work.  The workflow creates multiple files with various prefixes in tif format.  At the end of it all, the python script below does checks for bugs in some files, applies compression to some tif files, moves the needed files into a designated folder for training a CNN and then deletes files that are no longer needed.  If I run this in Spyder (and quit QGIS) it works fine.  But for the sake of simplicity and productivity, I'm trying to run it in the QGIS python console.  I've successfully installed skimage in the QGIS environment, so no problem with modules.  But, even if I remove all layers from QGIS, I still get permission errors that are presumably from a file lock imposed by QGIS. I need to reboot QGIS which cancels out any time gained when compared to running this in Spyder.   Is there a quick pyqgis command I could issue at the end of the code below to just release all locks? note that I don't even know a priori which layers are locked.

import os
import shutil
import numpy as np
from skimage.external.tifffile import imsave
from skimage import io

#############################################################
#User data input. Fill in the info below before running
PreProcessPath = "D:\\DeepRiverscapes\\PreProcess\\"
TrainPath = "D:\\DeepRiverscapes\\Predict\\"
RiverName = "StMarg"

################################################################
################################################################
#Reconvert the segments to a 1 band set of 1 label objects with integer labels without gaps

def RebuildSegments(Segments):
    FloatValues = np.unique(Segments)
    IntSegments = np.zeros(Segments.shape)
    S=1
    for n in range(0,len(FloatValues)):
        IntSegments[Segments == FloatValues[n]] = S
        S+=1
    return np.int32(IntSegments)

#Check for anomalies in the classification
def CheckClass(Class, name):
        if np.max(Class) > 7:
            print('WARNING some training class labels are greater than expected')
            print(name)
        elif np.min(Class) < 0:
            print('WARNING some training class labels are negative')
            print(name)
        elif np.count_nonzero(Class != np.int32(Class)):
            print('WARNING some training class labels are non-integers')
            print(name)

for i in range(0,99999): #goes to 99999 but do not use that many images!  It will just scan through the possible names
       ImPath = PreProcessPath + RiverName + format(i,'05d') +'_modified' + '.tif'
       ImPathIn = PreProcessPath + RiverName + format(i,'05d') + '.jpg'
       ClassPath = PreProcessPath + 'R_VECT_' + RiverName + format(i,'05d') +'_modified'+'.tif'
       SegPath = PreProcessPath + 'SEG_' + RiverName + format(i,'05d') +'_modified' + '.tif'
       ImPathOut = TrainPath + RiverName + format(i,'05d')  + '.jpg'
       ClassPathOut = TrainPath + 'CLS_' + RiverName + format(i,'05d') +'.tif'
       SegPathOut = TrainPath + 'SEG_' + RiverName + format(i,'05d') + '.tif'
       if os.path.exists(ClassPath):

           Class = io.imread(ClassPath)
           CheckClass(Class,ClassPath)
           imsave(ClassPathOut, Class, compress=7)
           Segments = io.imread(SegPath, as_grey = True)
           Segments_int = RebuildSegments(Segments)
           imsave(SegPathOut, Segments_int, compress=7)
           #Cleanup and move operations
           shutil.move(ImPathIn,ImPathOut)
           os.remove(SegPath)
           os.remove(ImPath)
           os.remove(ClassPath)


texts: Are "gridded data" and "raster data" the same thing?

texts:
I'm working on a project that will involve both "raster data" (aerial imagery and DEMs) and "gridded data" which I understand to be other kinds of numerical measurements taken at grid points, and which could be completely non-image related (eg, population density).

I'm having trouble finding much on "gridded data", so I'm wondering if it's just a different term for raster data? Is there any fundamental distinction? Are there different tools (and particularly servers) for these two types?

(If relevant, this is in the context of Australian state governments...)


texts: Filled Polygons using ggplot in R not working

texts:
The following code works well for drawing polygons

SP <- SpatialPolygons(spd)
centroids  <- coordinates(SP)
x  <- centroids[,1]
y  <- centroids[,2]
SP1  <- SpatialPolygonsDataFrame(SP, data = data.frame(x = x, y = y, z = z, row.names = row.names(SP)))
ggplot(data = SP1, aes(x = long, y = lat, group = group)) + geom_path()

Instead if I change the last two lines of the code to fill these polygons based on some value vector it does not work:

penetration <- runif(n = length(SP), min = 1, max = 10)
SP1  <- SpatialPolygonsDataFrame(SP, data = data.frame(x = x, y = y, z = z,  value = penetration, row.names = row.names(SP)))
ggplot(data = SP1, aes(x = long, y = lat, group = group)) + geom_polygon(aes(fill=SP1$value))

spd is created with data from the following link:

data source: https://www.google.com/maps/d/viewer?mid=1mf9TQY3bocarO4SaIXAtco_C1LE&hl=en_US

I then extracted coordinates from the kml file... which is then into a list of Polygons so that I could draw it using some plotting function...

tkml <- getKMLcoordinates(kmlfile="Pune all Electoral Wards.kml", ignoreAltitude=T)
len <-length(tkml)

del <- NA
for(i in 1:len)
{
  if(length(tkml[[i]]) == 2)
    del <- c(del, i)
}
del <- del[-1]
tkmlCleaned <- tkml[!tkml %in% tkml[del]]
size <- length(tkmlCleaned)
sr <- vector(mode = "list", length = size)
srs <- vector(mode = "list", length = size)
spd <- list()
for(j in 1:length(sr))
{

  sr[[j]] <- Polygon(tkmlCleaned[[j]])
  srs[[j]] <- Polygons(list(sr[[j]]), j)
  spd <- c(srs[[j]], spd)
}


texts: Nearest Neighbor distance frequency distribution

texts:
I have a dataset of 6455 points. Is there a way to produce the frequency distribution of nearest neighbour distances in the data set in ArcGIS 10.3?
The Average Nearest Neighbor tool only produces the mean distance, not the whole distribution.


texts: Missing Symbology Categories after performing join using ArcObjects?

texts:
I am trying to join a table to a feature class in , and use the joined fields to create custom symobology based on the joined fields as the picture 1 below  for Arcmap using ArcObjects. I am successfully joining the tables, and I can see the new fields in the layer when I open it's attribute table, however once successfully joined, and the code runs to specify the rederer I check the symbology and the entire category criteria gets hidden as picture 2. What could be causing this ? It is very strange, when I remove the code where the layer and table get joined, the symbology becomes correct using the 3 values I specified and it doesn't disappear.

[

Here my code

   protected override void OnClick()
    {

    // code that gets all the proper objects
    //..
       JoinTables(table,ftrClass,ftrLayer);
       var geoLayer=(IGeoFeatureLayer)ftrLayer;
       CreateRenderer(geoLayer);

    }

      // method to join table and ftr class layer
     public void JoinTables(  IStandaloneTable pStTable,IFeatureLayer ftrLayer, IFeatureClass ftrClass){
       var pTabCollection = (IStandaloneTableCollection)arcMapApp.Document.FocusMap;
        if (pTabCollection.StandaloneTableCount == 0) {
          return;
        }
        var pDispTable2 = (IDisplayTable)pStTable;
        var pTTable = pDispTable2.DisplayTable;
        IMemoryRelationshipClassFactory pMemRelFact;
        IRelationshipClass pRelClass;
        pMemRelFact = new MemoryRelationshipClassFactoryClass();
        pRelClass = pMemRelFact.Open("JoinedLayerAndTable", (IObjectClass)pTTable, "CL_ID",
          (IObjectClass)ftrClass, "ID", "forward", "backward", esriRelCardinality.esriRelCardinalityOneToMany);
        var dispRC = (IDisplayRelationshipClass)ftrLAyer;
        dispRC.DisplayRelationshipClass(pRelClass, esriJoinType.esriLeftOuterJoin);
  }

   public void CreateRenderer( IGeoFeatureLayer ftrLAyer){

        IUniqueValueRenderer renderer = new UniqueValueRendererClass();

        renderer.FieldCount = 3;
        renderer.Field[0] = "CATEGORY_NAME";
        renderer.Field[1] = "TYPE_NAME";
        renderer.Field[2] = "UUID";
        renderer.FieldDelimiter = "|";
        renderer.ColorScheme = "Custom";
        renderer.FieldType[0] = true;
        renderer.FieldType[1] = true;
        renderer.FieldType[2] = true;

        ILineSymbol smb = new SimpleLineSymbolClass();
        smb.Width = 2;
        IRgbColor clr = new RgbColorClass();
        clr.Red = 255;
        clr.Green = 0;
        clr.Blue = 0;
        smb.Color = clr;

        renderer.AddValue("XX1|INSUR|DWA", "www", (ISymbol)smb);
        pGeoFeatureLayer.Renderer = (IFeatureRenderer)renderer;
        var pUID = new UIDClass();
        pUID.Value = "{683C994E-A17B-11D1-8816-080009EC732A}";
        pGeoFeatureLayer.RendererPropertyPageClassID = pUID;
        arcMapApp.Document.UpdateContents();
        arcMapApp.Document.ActivatedView.Refresh();
  }


texts: Calculating distance between point layer to n nearest polylines using QGIS?

texts:
In QGIS I am trying to calculate the distance from a point layer to a polyline layer for the n nearest polylines. I have figured out how to calculate this distance to the centroid of the polylines using the point to hub tool, but I need the distance to the edge of the closest polylines - not centroids. I can calculate the nearest polyline by dumping the polylines to nodes and calculating to nearest point, but that method falls apart for the n nearest calculation.

Any thoughts?


texts:
I am new to QGIS and want to import and view the Swedish .rik.

I am looking for simple instructions how to import the .rik files into QGIS and then view and print the maps.


texts: Adding Layer using esriRequest Response in ArcGIS API for JavaScript?

texts:
For various reasons I need to add layers to a map object using a response returned by an esriRequest to an ArcGIS Server. One reason being the data will change constantly and hard coding FeatureLayer URLs would prove inefficient.

I have successfully been able to return a JSON object containing attribution that appears to be the same data one would see on the layers REST node.

 maps = [];

require([
    "esri/map", "esri/layers/FeatureLayer",
    "esri/dijit/Print", "esri/tasks/PrintTemplate", "esri/request",
    "esri/config","esri/dijit/BasemapToggle",

    "dojo/_base/array", "dojo/dom", "dojo/parser",  "dojo/domReady!"

], function (Map, FeatureLayer,Print, PrintTemplate, esriRequest,      esriConfig,BasemapToggle,arrayUtils, dom, parser){

parser.parse();

console.log("Map Start");

esriConfig.defaults.io.proxyUrl = "http://........./proxy.ashx";

maps.addRespLyrs = function addRespLyrs (layer) {
    if (layer.type == "Feature layer") {
        console.log(layer.name);
        maps.map.addLayer(layer);
    }
};

maps.handleMapInfo = function handleMapInfo (resp) {
    console.log("handleMapInfo called");
    maps.layers = resp.layers;
    console.log(maps.layers);
    for (i in maps.layers) {
        var layer = maps.layers[i];
        maps.addRespLyrs(layer);
    }
};

maps.layerRequest = function layerRequest () {
    var mapInfo = esriRequest({
        "url": "https://........./Application Name/MapServer/layers",
        "content" : {"f":"json"},
        "handleAs" : "json",
        "callbackParamName" : "callback"
    });
    mapInfo.then(maps.handleMapInfo);
};

maps.map = new Map("cpCentCent", {
    basemap:"topo",
    center: [-98.5795, 39.8282],
    zoom: 5
});

console.log("Map End");

When attempting to add the data using the addLayer() method I recieve an error thrown from deep with in esri's API.

  init.js:199 TypeError: a.id.match is not a function(…) "TypeError: a.id.match is not a function

Any thoughts as to what could be going wrong here or how to create a layer from my JSON?


texts: Zoom to a plot with GeoPandas based on data from CSV and shapefile

texts:
Here is my code:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import descartes
import geopandas as gpd
from sklearn.cluster import KMeans
from sklearn.preprocessing import MinMaxScaler
from sklearn.decomposition import PCA
from shapely.geometry import Point, Polygon
from geopandas import GeoDataFrame as gdf
%matplotlib inline

street_map = gpd.read_file('C:\\Users\\user\\Downloads\\raccoon_island\\raccoon_island.shp')

fig,ax = plt.subplots(figsize = (15, 15))
street_map.plot(ax=ax)

df = pd.read_csv('F:\\Anaconda\\docs\\OneDrive_1_9-29-2020\\11.csv')
df = df[['longitude_start','latitude_start']]
crs = {'init': 'epsg:4326'}
df.head()

geometry = [Point(xy) for xy in zip( df["longitude_start"], df["latitude_start"])]
geometry[:3]

geo_df = gpd.GeoDataFrame(df, crs = crs, geometry = geometry)
geo_df.head()

fig,ax = plt.subplots(figsize = (15, 15))
street_map.plot(ax = ax, alpha = 0.4, color='grey')
geo_df.plot(ax=ax, color='red', markersize=20, marker = '+', label= 'pos')
plt.legend(prop={'size': 15})

I have tried the following but it only seems to make the outer box change size. I would like to zoom in on the result instead.
fig,ax = plt.subplots(figsize = (30, 30))
street_map.plot(ax = ax, alpha = 0.4, color='grey')
geo_df.plot(ax=ax, color='red', markersize=20, marker = '+', label= 'pos')
ax.set_ylim(auto=True)
plt.legend(prop={'size': 15})

ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax))

minx, miny, maxx, maxy = gdf.total_bounds
ax.set_xlim(minx, maxx)
ax.set_ylim(miny, maxy)


texts: Trying to modify records from one shapefile to another gives Error 999999?

texts:
I am trying to update the atrribute data of a feature class with respect to another one by matching a common field. But it is returning an error for the setValue command. My code looks like:

import arcpy
fc="D:/Geonames_CA.gdb/Export_Output"
fs="D:/abc.gdb/Export"
rows = arcpy.SearchCursor(fc)
fldList=arcpy.ListFields(fc)
rws = arcpy.SearchCursor(fs)
flList=arcpy.ListFields(fs)
for row in rows:
    for fld in fldList:
        a = row.getValue("PLACE_NAME")
        a_xmin=row.getValue("XMIN")
        a_xmax=row.getValue("XMAX")
        a_ymin=row.getValue("YMIN")
        a_ymax=row.getValue("YMAX")
for rw in rws:
    for fl in flList:
        b = row.getValue("PLACE_NAME")
        b_xmin=row.getValue("XMIN")
        b_xmax=row.getValue("XMAX")
        b_ymin=row.getValue("YMIN")
        b_ymax=row.getValue("YMAX")
        if b==a:
            cursor = arcpy.UpdateCursor(fs)
            for r in cursor:
                r.setValue(b_xmin, a_xmin)
                r.setValue(b_xmax, a_xmax)
                r.setValue(b_ymin, a_ymin)
                r.setValue(b_ymax, a_ymax)
                cursor.updateRow(r)
            del r, cursor

I am only getting this error:

  Traceback (most recent call last):   File
  "C:\Python27\ArcGIS10.1\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py",
  line 326, in RunScript
      exec codeObject in main.dict   File "D:\Script2.py", line 25, in
      r.setValue(b_xmin, a_xmin)   File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line
  1007, in setValue
      return convertArcObjectToPythonObject(self._arc_object.SetValue(*gp_fixargs(args)))
  RuntimeError: ERROR 999999: Error executing function.


texts: FME - process only matching name datasources

texts:
I'm trying to find out how to deal with the below scenario with FME 2012 on board.
I have a bunch of geotiff's tiles (A) and a bunch of others geotiff's tiles (B) that I would like to mosaick together.
The name of each tile from (A) match eaxactly one tile in (B).
I would like the process to be run tile by tile due to the amount of the tile I have to process (roughly 13k).

What whould be the scenario with FME to:
- take a pair of tiles with exact names from A and B,
- run a workbench,
- output the result,
- take another pair of tiles with exact names from A and B,
...
and so on.

The problem is that I actually don't know how to force FME to read a tile A then find a coresponding tile in B (not reading all tiles from B) and run a workbench for the matching pair.

Thanks in advance for help with that.
G.


texts: SelectFeature and ModifyFeature in OpenLayers 2.13

texts:
After lot of experimenting, I have to give up and ask: has anyone had success in using both SelectFeature and ModifyFeature in OL 2.13 (NB, there were changes to ModifyFeature in 2.13, compared to 2.12)?

What I want to do: ModifyFeature should work on it's regular way (user can drag points; and click on lines will bring "handles" visible). In addition, I want to catch select and hover events.


texts: Contours with negative values using the model example in QGIS?

texts:
I have a data set of water level ranging from -6 m below sea level to +30 meters.
But when I use the model example for creating contours from points, the minimum contour value is zero (0).
What can be the reason for that?

(I'm using QGIS 2.4.0)

UPDATE:
I've tried to run the model in QGIS 2.6 and get .
It works fine when running another data set with only possitive values.
What is the meaning of the error?


texts: Get the different values assigned to a zip code

texts:
I have a map to which I added a zip code layer and another layer which contains values assigned to every square kilometer within each zip code. I'm trying to find out how many different values there are in each zip code and what they are.

So for example, if a zip code has 10 values assigned to it in total, 5 of them being 1.5, 3 of them being 2.1 and the remaining two being 6.4, this zip code basically has 3 different values, namely 1.5, 2.1 and 6.4. That is what I want to know. Is there a way to do this without manually checking all the values in every zip code?

This is a screenshot of my data. The value that I'm working with is under the column BI. I need a table which gives me a list of the long/lat in every zip code and the corresponding values. I can't do anything on excel yet since I don't know which lat/long goes with which zip code.


texts:
I am trying to compute population affected by various distances, during earthquake events. The earthquake point events update every minute and I need to find a way how to create automatic buffers and get real time statistics. At the moment, I have point data and have used geometry generator polygons to create the buffer affect for 5km, 10km and 20km (as can be seen via attachment). What is the best way to calculate zonal statistics using these polygons created through geometry generator?
I am also using Print Composer Atlas to export these maps in real time and would like to add the affected population for each buffer through a legend key on the map below. Is there any way I can do this?


texts: How do I configure JVM to run in server mode when using Geoserver as a Windows service?

texts:
I have been reading this article about putting Geoserver into a production environment.

Section "4.1 Java Virtual Machine" is confusing me a bit. This section mentions that in a production environment Geoserver should be using the JVM in "server mode" and it also shows you how to check what mode your JVM is currently in. As it turns out mine is running in "client mode", which appears to be default for windows.

I thought that I could get Geoserver to use the JVM in "server mode" by adding this line to my wrapper.conf file.

wrapper.java.additional.4=-server

After I added that line, I restarted my service (even bounced the server on the second attempt) and checked the JVM again by going "java -version" but it still appears to be in "client mode".

I know my JVM can perform server mode because this "java -server -version" returns results like in the article.

Have I done the right thing? Or have I missed a step somewhere?

Windows 2003 - Geoserver 2.1.4 - Java version 1.6.0_23

Thanks

Ando


texts: Importing data into a project with a custom CRS incorrectly projected

texts:
I created a custom CRS for a project I am working on.  I then started a new project, set to my custom CRS and turn on "on the fly" transformation.  I import Rail data that I know is in my custom CRS, and it imports as expected.  However when I import Road data from our library that is in a different CRS, this data goes to a totally different area.

I tried the reverse situation.  I opened the Road data in a new project, set the CRS to this layer and turn on "on the fly" transformation.  When I import the Rail data (which is in my custom CRS) it imports fine.

Is it possible to create a project based on a Custom CRS?

Ultimately I want to import the data I am interested in, convert it to my custom CRS so I can test out some of the grass functions (which I understand needs to be in the same CRS), and use qconsolidate to store it is a separate location.


texts: Convert a RGB 3-band Tiff to a singleband colored image in ArcGIS

texts:
I have a RGB 3-band GeoTiff which gives me a nice looking when loaded in ArcGIS. Now i want to convert it to a simple colored image. My goal is to have an image which can be displayed the same in different software. (I know that i can change symbology of raster in every software, but i don't want to do that.) I have no analysis or processing need for this image, and only its view is important for me. How can i do that in ArcGIS?
What i have done: using Copy Raster tool, but unfortunately 'RGB to colormap' option is grayed out for this raster. I also tried Raster Calculator to add band values together, but output is not colored.


texts: No method or default for coercing "sf" to "SpatialPoints"

texts:
I keep getting the below error message when I try to use R in QGIS.  For example Monte Carlo or Quadrant.
I use a point shape file with the x and y coordinates as centroids within the attribute table.
For example
##Point pattern analysis=group
##Layer=vector point
##Simulations=number 100
##Optional_plot_name=string
##showplots
library(spatstat)
library(maptools)
sp <- as(Layer, "SpatialPoints")
sp <- as(sp, "ppp")
e <- envelope(sp, Kest, nsim = Simulations)
>e
plot(e, main = Optional_plot_name)

Error Message
Error in as(Layer, "SpatialPoints") :
no method or default for coercing "sf" to "SpatialPoints"
Calls: as -> .class1 -> as

Info on the point shapefile as follows
Name
LandScanJoinCTAT
Path
G:\Article\GIS\LandScan\Master\LandScanJoinCTAT.shp
Storage
ESRI Shapefile
Comment

Encoding
System
Geometry
Point (Point)
CRS
EPSG:4326 - WGS 84 - Geographic
Extent
17.0590112291777913,-34.4981314092633582 : 32.5317232060707155,-22.4313408665917606
Unit
degrees
Feature count
234


texts:
I have a series of data points that are placed out along an arbitrary grid and I have the known UTM of the datum for this grid.

I'm sure there is a way to snap all the arbitrary points so they are placed out in real UTMs based on the known UTM of the datum, but because it's early in the morning and I haven't had my coffee yet, I can't seem to figure it out.


texts:
I'm using Landsat 8 imagery, and the Indicar Landsat Processing Tools to produce NDVI imagery from the Landsat 8 source. That works fine.

I then have a bounding box expressed as lat/lon coordinates, and I want to crop the NDVI image using gdalwarp:

gdalwarp -s_srs "+proj=utm +zone=32N +datum=WGS84" -t_srs WGS84 -te 8.45586 56.55605 8.62229 56.78750 LC8_197-021_20130502_LGN01_ndvi.tif cropped112.tif

The original (left) vs. result (right) looks like this:

I'm sure it has to do with the value of t_srs. I tried not to specify neither s_srs or t_srs which fails with: ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.

Probably because it expects bounding box in something different that lat/lons? Meters?

gdalinfo output for the original NDVI image:

Driver: GTiff/GeoTIFF
Files: LC8_197-021_20130502_LGN01_ndvi.tif
       LC8_197-021_20130502_LGN01_ndvi.tif.aux.xml
Size is 7921, 7711
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (364785.000000000000000,6311415.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  364785.000, 6311415.000) (  6d46'41.90"E, 56d55'35.84"N)
Lower Left  (  364785.000, 6080085.000) (  6d53'38.18"E, 54d50'58.48"N)
Upper Right (  602415.000, 6311415.000) ( 10d40'59.05"E, 56d56' 6.12"N)
Lower Right (  602415.000, 6080085.000) ( 10d35'43.57"E, 54d51'26.50"N)
Center      (  483600.000, 6195750.000) (  8d44'15.63"E, 55d54'24.87"N)
Band 1 Block=7921x1 Type=Float32, ColorInterp=Gray
  Min=-0.232 Max=0.893
  Minimum=-0.232, Maximum=0.893, Mean=0.138, StdDev=0.210
  Metadata:
    STATISTICS_MAXIMUM=0.89341741800308
    STATISTICS_MEAN=0.1375360206514
    STATISTICS_MINIMUM=-0.23159252107143
    STATISTICS_STDDEV=0.2095535427118

gdalinfo output of the warped image:

Driver: GTiff/GeoTIFF
Files: cropped112.tif
Size is 424, 589
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (8.455859999999999,56.787500000000001)
Pixel Size = (0.000392523584906,-0.000392954159593)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   8.4558600,  56.7875000) (  8d27'21.10"E, 56d47'15.00"N)
Lower Left  (   8.4558600,  56.5560500) (  8d27'21.10"E, 56d33'21.78"N)
Upper Right (   8.6222900,  56.7875000) (  8d37'20.24"E, 56d47'15.00"N)
Lower Right (   8.6222900,  56.5560500) (  8d37'20.24"E, 56d33'21.78"N)
Center      (   8.5390750,  56.6717750) (  8d32'20.67"E, 56d40'18.39"N)
Band 1 Block=424x4 Type=Float32, ColorInterp=Gray


texts:
I am using ER Mapper. But, I'm still confused how to detect and mask the clouds.

I want to reduce or eliminate the cloud to perform supervised and unsupervised classification.


texts: Take from multilinestring the start and end points

texts:
I have to use the A* algorithm and my geometry is in multilinestring and I want to take stratpoint(x1,y1) and endpoint(x2,y2) .
I use this queries but the columns after the execute is still empty

UPDATE roads SET x1 = ST_x(ST_startpoint(geom));
UPDATE roads SET y1 = ST_y(ST_startpoint(geom));
UPDATE roads SET x2 = ST_x(ST_endpoint(geom));
UPDATE roads SET y2 = ST_y(ST_endpoint(geom));

after the execution give this message

  Query returned successfully: 2584 rows affected, 101 ms execution
  time.

Then I make a select and the columns are empty

I searched the net and I didn't find anything, I don't know what I can do.


texts: Focal / Moving Window statistics - open source

texts:
I am currently running a process that involves the creation of many "Focal Statistics" rasters in ArcGIS. These rasters use varying-radius circular search areas and calculate three different statistics: Mean, Sum, and Standard Deviation (depending on the input data source). They also must be created in a programmatic manner - either from Python, R, or Bash/Command Line.

I would like to move this to an open source alternative but am having trouble finding an open source equivalent to ArcGIS's "Focal Statistics". The raster::focal() function in R is terribly slow, and for the size of my data it will not work. SAGA GIS has the 'simple filter' tool, but that will only calculate the Mean, not the other two statistics.

Does anyone have an alternative that replicates the functionality of "Focal Statistics" in a reasonably fast, free, and programmatic manner with a circular search window?


texts: I'm looking for a DEM of Russia projected in WGS84. Anyone have any ideas on where to look?


texts:
I have seen that this issue has already been reported in the past but for previous QGIS versions. So, I decided to raise it again and, in the case it would be considered a repetition, I excuse.
I am using QGIS 3.16.0-Hannover and QuickMapServices (QMS) version 0.19.11.1.
QMS was working while, all of a sudden, it stopped showing basemaps. I mean every basemaps available for QMS.
Maybe it's my fault because I could have changed some settings but I really do not know what to do.


texts: Loading Natural Earth *.prj file into QGIS 2.10.1?

texts:
I have found solutions for this problem in other versions of QGIS, but none that apply to the newest version 2.10.1. Whenever I try to import a Natural Earth Raster projection that is a prj file I get an Invalid Data Type error.

I have tried using the old solutions to specify the file I'm importing, but its no longer supported.

Im trying to load the Natural Earth data set "NE1_HR_LC_SR_W.prj" @ http://www.naturalearthdata.com/downloads/10m-raster-data/10m-natural-earth-1/


texts: Problems updating geometry of multipart and donut polygons with Arcpy

texts:
I am trying to update the geometry of some features (polygons) to allow stretching / shrinking in the X and Y planes, but am struggling to get this to work for multipart polygons and for those with an inner ring (donut polygons). I am struggling to find a solution online and with weak Python skills and suspect something simple is wrong in my code - possibly indentation in the final part of the program?

Two issues arise in this code- 1.) multipart polygons are only updated for their first part and therefore become single part, and 2.) the inner ring is lost in the donut polygon.

Attaching a before processing and after image to show the issue, shapefile and code. Looking at the after image (cyan polygons), my blue (multipart) polygons are not replicated properly, and my donut polygon (red) becomes a single area with no inner gap.

Is this just indentation, or do I need more looping to deal with the multipart & inner ring logic correctly?

Code:

import arcpy

def stretch_features(in_features, x_stretch=None, y_stretch=None):

    with arcpy.da.UpdateCursor(in_features, ['OID@', 'SHAPE@XY', 'SHAPE@']) as cursor:
        for row in cursor:

            Xc = row[1][0]
            Yc = row[1][1]

            print("Feature {}, at {}, {}".format(row[0] ,Xc, Yc))

            geometry = row[2]

            partnum = 0

            for part in geometry:

                print("Part {}:".format(partnum))
                pts = geometry.getPart(partnum)

                pt_array = arcpy.Array()

                for pnt in pts:
                    if pnt:
                        #print("old {}, {}, new {}, {}".format(pnt.X, pnt.Y, Xc-((Xc-pnt.X) * x_stretch or 1), Yc-((Yc-pnt.Y) * y_stretch or 1)))
                        pnt.X = Xc-((Xc-pnt.X) * x_stretch or 1)
                        pnt.Y = Yc-((Yc-pnt.Y) * y_stretch or 1)
                        pt_array.add(pnt)
                    else:
                        print("interior ring")
                        print("old inner {}".format(pnt))
                        pt_array.add(pnt)
                        #pt_array.add(None)

                partnum += 1

            new_polygon_geometry = arcpy.Polygon(pt_array)

            row[2] = new_polygon_geometry

            cursor.updateRow(row)
    return

in_features = r"D:\ReShapeFeatures\TestPolygons.shp"

stretch_features(in_features, x_stretch=.75, y_stretch=.75)


texts: Change unit in Split tool on ArcGIS Pro?

texts:
I'm working on a project with a CRS in meters, map units and display units are also in meters, but Split tool (from Modify features) is in feet:

Is there a way to change unit to split based in a distance in meters?


texts:
I am using the scale bar of the QGIS composer (ver 1.8) and I don't understand how the whole thing works.

Where can I set the units properly? The layers have idential CRS (WGS84), the settings for map tools is in meters.

I have tried to get it correctly by changing segment size and map units per bar unit, but it doesn't work. However, I have to admit that I don't understand what these two fields actually mean.