texts: Dissolve polygons in Single-parts and calculating statistics for those single-part polygons only

texts:
I have a polygon shapefile. Which have some numeric attributes (A, B). I clustered the shapefile in five groups (labels) based on an attribute(A) and want to dissolve them while aggregating other attributes (A, B).
Code for that is:
import geopandas as gpd
# suppose clustered_df is a gdf obtained after clustering
db_clustered_dissolve = clustered_df .dissolve(by='labels'
            , aggfunc= {'labels':'count','A': 'sum', 'B': 'sum'})

The above code dissolves the dataframe by labels column and create Multi-part polygons.
I know that if I want only single-part polygons then I can use explode() on the dataframe db_clustered_dissolve but I can not accuratly calculate A, B and count for those resulted single part dissolved polygons because every single-part polygon will have same values for A, B and count as their corresponding Multi-part polygons.
I'm having this idea that I can use something that only single-part polygon would be created while dissolving and statistics for those could be calculated on the spot.


texts: Address Locator Only Fails with Suggestions Enabled

texts:
I am currently trying to create a General Single Field Address Locator in ArcMap 10.3.1. I can create the locator fine without the suggestion option enabled, but with it’s on it gives me the Error 000042. I have already removed special characters, and any duplicate names. I also tried saving off our network, in case of permissions issues to no avail. The error occurs when it is running the “Building suggestions index (merging data)…”

Are there any ideas out there what might be causing this?


texts: Realtime imagery from Google Earth Engine

texts: Can we use Google Earth Engine to get realtime (not later than a week is also fine) satellite images?


texts: How to find and apply reprojection parameters?

texts:
I have a data set in one projection, and a shapefile in another, and I am trying to combine the two. I am using Python. I apologize for posting a "help me transform coordinates" question on this site - there are so many.

I assume the way to do this conversion is with proj4. I am unable to find/understand the proj4 string needed to do this conversion.

The shapefile.prj contains this:

PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",90.0],UNIT["Meter",1.0]]

Is there a simple way to find the proj4 code for this?

And I have a netCDF file that I can plot correctly with the following Python code:

map = Basemap(llcrnrlon=lon[0,0],
              llcrnrlat=lat[0,0],
              urcrnrlon=lon[560,300],
              urcrnrlat=lat[560,300],
              #rsphere=(6378137.00,6356752.3142),\
              resolution='l',
              area_thresh=1E6,
              projection='lcc',
              #lat_1=90, lon_0=0)
              lat_1=71.0,
              lon_0=-39.0)

I would like to keep the Basemap (python) projection, and convert the shapefile contents to that projection.

Edit

Further research and I found this post: Shapefile PRJ to PostGIS SRID lookup table? with a helpful bit of code. Running that and I find the following:

Proj4 is: +proj=stere +lat_0=90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
EPSG is: None

Edit

I found out I can determine the basemap proj string with:

map.proj4string

I'm now using the advice here http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/ to try to convert from one to the other...


texts: Image does not display after clipping raster using ArcGIS Desktop?

texts: I am using ArcGIS 10. I am trying to clip an NLCD layer from a vector polygon. I have used the Clip Raster tool in Data Management and Extract by Mask tool in Spatial Analyst among other things. My problem is that the processing completes without an error message but the image does not appear. I have zoomed to the layer and I have also checked in ArcCatalog.


texts:
I am new to GIS.

Which program do you recommend I use to view, mark, and then print .rik map files in Linux?


texts: Can I use geoserver to bring sql layer as vector so I can edit it?

texts:
I've updated geoserver with SQL Store and published a table of geometries.
When using WMS I get the LineStrings as img so I can't edit them.
I've read that WFS bring the data as vector (features) but I wasn't able to find example how to do it (When I preview layer in OpenLayers it  comes as WMS only)

How can I do it ?


texts:
I'm trying to reproject this dataset from Austria Statistik Geodata.

Data comes in "MGI_Austria_Lambert", and I'm trying to reproject to WGS84

What I'm doing is Saving layer as... and selecting EPSG:4326, but the layer is still in the same place.

The other option i'm trying is changing the original .prj file with other one from a layer in EPSG:4326, but I have the same result.

I don't what can I do to solve it.

I've also tried setting the projection with:Set SRC to layer, and selecting 4326, but no changes.


texts:
I’m working with a fresh install of QGIS 3.0.3 on Ubuntu 18.04 and am not able to get my XYZ Tiles to render. When I add an XYZ Tile to the Map View it doesn’t appear. When I look at the XYZ Tile properties they look normal, as if the coverage was part of the project. If the XYZ tile is the first coverage that I open in a new project the project projection is set to the tile’s projection.

Although it’s quite possible that it’s a user problem, I followed several tutorials and the process seems quite simple. Any suggestions on ways to fix this or troubleshoot are welcome.


texts: How to edit specific geometry on leaflet layer

texts:
I am starting working on a new, quite complicated functionality and I want to ask for some guidance regarding the way to implement it.
The final aim is to give to the user the possibility to edit data of a WFS layer.
Right now I have completed the editing and updating of the non spatial data.

More specifically, I used DataTables to present all the attributes/values of the WFS layer (got the values through OWSLib). I then updated the values by sending a POST request (using python Requests).

That was quite straightforward. Now I want to add geometry editing functionality.

As its shown in the screenshot, I want to add a Edit Geometry button next to Edit button.
The idea is that on button click, I get a modal with the map (I will use Leaflet) and all the layer loaded.

I see there lots of Leaflet plugins for drawing. But what I want is to be able to edit only the feature of the specific row (e.g. with fid = 3).
So this is my challenge. How to show the whole layer but giving editing permissions only to the feature of the selected row?

Is out there one of the leaflet plugins supporting this? Or I should customize it?
Would you suggest a specific leaflet plugin for drawing/editing geometries?


texts: I have selected certain values in a raster (roughly 1000 cells), and would like to export them to be able to create a polygon. Converting it to points first then to polygons does not work. Can anyone suggest a better method?


texts:
I am a R user and quite new in GIS.

From Adding Attribute data to shapefile? I know how to add/edit attribute data (say population) to shapefile using R.

Now, I have some shapefiles with attribute data, and I know how to make a static map, but I would like to know if there is a way to make an interactive map in R, e.g. allow zooming, hover a region will display the corresponding attribute data in R.

If not, I have heard of a free and open-source GIS software QGIS, is it possible to do it using that instead?


texts: Export topology tables to txt files using ArcObjects, generating directed graph, and finding intersections

texts:
I want to export topology tables to  text files:

for points:
1  id  x   y

for lines:
line id   fompoint.x  from point.y  topoint.x  topoint.y

for polygons:
id   line id 1, line id 2...

2   arc list
3   node list

Using ArcObjects VBA.


texts: Huge memory consumption rendering geojson files

texts:
I have a road network graph saved in geojson - size about 150MB. When I load it to QGIS as the only layer, QGIS (which uses about 140MB at this point) starts allocating memory up to about 8GB. Another memory (GBs) is then allocated when zooming in/out or moving the view quickly hitting all available RAM (about 12GB). Is this normal QGIS behavior or is it a bug/memory leak?

EDIT: tested on Windows on the same data - memory usage does not exceeded 2GB and decreased after full load of the layer. So it looks like the problem is linux/ubuntu specific.


texts:
I'm a fairly new user to QGIS and was trying to generate a map of Iceland. However my grid scales are not matching up. The base map I'm using is Esri National Geographic which uses CRS- EPSG:3857 - WGS 84 / Pseudo-Mercator - Projected

The rough width of Iceland is 500km but according to the grid boxes QGIS generated in Print Layout its about 1000km. In grid properties they are set to be 5000000m x 500000m. And the scale bar is not accurate either.

Screenshots attached.


texts: How to resolve ArcGIS for Desktop Error (SDE has lost its connection with underlying DBMS[ORA-03114:not connected ORACLE])?

texts:
How to resolve this issue when user leaves the system idle for one hour and again trying to connect Database, It encounters the Connection error message and then user should close the ArcGIS and restart it again.
Is there any temporary solution to make the DBMS connection active for one user if needed to stay more than one hour idle and work again?
Thanks


texts: Copy TIGER geocoder data to another machine

texts: I have downloaded census data and I want to copy to a different machine, which files should I move and how can I create the tables from this? since I will not be using the regular loader scripts.


texts: Merge feature classes in feature datasets using ArcPy?

texts:
I'm trying to merge a group of like feature classes which are contained within separate feature datasets. I can merge feature classes within a file geodatabase without having issues by creating a list of the feature classes and using Merge_management. The trouble comes when I need to do the same for feature classes contained in feature datasets.

My code so far, looks like:

Output = filepath for new fc
for datasetList1 in arcpy.ListDatasets("","Feature")+ [""]
    PointsFC in arcpy.ListFeatureClasses("*_PointName*","Point",datasetList1):
    for points in PointsFC:
    arcpy.Merge_management(points,Output)
        print "Merged Points"

The error message I get states:

  PtsMerge already exists.
  Failed to execute (Merge).


texts: I have a polygon layer that consists of many different polygons; some of them are neighbouring and some of them are not. How could I merge the neighbouring polygons?


texts: Joining multible rows in csv file to a single row in shp file using ArcGIS

texts:
I have two tables (one is actually a shapefile and the other one is a comma separated file).

I really need your help with joining the csv file table (which has many rows for one ID number) with the shapefile attribute table, so when I'll identify something in the map, it'll show me all those rows.

I tried joining, I also tried importing these to geodatabase and I tried to change the field type of the csv attribute table (because it didn't let me to join otherwise). But it didn't work, those other two rows went still missing.

My steps doing this with geodatabase files.

First I made a geodatabase:

Then I imported the shapefile and the tabel:

I dragged both of these files to arcgis and joined the tables

And when I searched the ID I know has 3 rows in the original table, I saw that it hasn't added them here in the joined tables (here is only the first row

Related tables

There are none related tables after I've made the join.

But after I made the join and relate also, I could choose to see them in the table.

And I find a way to see the relates also in the indetify box, I had to open that "tree".


texts:
I've read how to use the QGIS Browser to add layers into the project from within zip files, but this doesn't seem to be available with processing dialogs.

I'm trying to merge two raster files using Raster>Miscellaneous>Merge. For my input files it only recognises unzipped files in the folder. If I choose the zip folders as inputs it gives and error that this is not a recognised file type. Is there a way of taking files into processing dialogs direct from zip folders?


texts: Specifying grid width and height dimensions with ArcMap custom toolbox parameter

texts:
I am using I am writing a script to create a tool (in ArcMap 10.7.1) that involves overlaying an agricultural field with a specific grid size. I need the grid to contain a custom width and height of the individual cells, as specified by the tool user. The input will look like this:

#User-defined input, stored as an integer
cellSize = arcpy.GetParameterAsText(0)

#Generate grid
polygonWidth = cellSize
polygonHeight = cellSize

The grid will be generated with the arcpy.GridIndexFeatures_cartography tool.

What I don't know how to do is to add the measurement units to the end of the tool parameter. I want the numeric entry that my users provide to be followed by "feet."

When I upload my script to the toolbox, is there an option that will allow me to specify the measurement units, so that I will be able to generate a proper grid?


texts: Renaming 1,200 rasters in File Geodatabase using ArcPy

texts:
I've been tasked with renaming a large amount of raster files within a file geodatabase (ArcGIS style).

I'm taking weather data which has a name format as such: "G"mmddhhhh
I'm trying to add a year tag in between the G and month: "G"yyyymmddhhhh

Here is the function that I've written:

# ----------------------------------------------------------------------------- #
@timeit
def renameRasters(FGDB):
    origdir = arcpy.env.workspace
    arcpy.env.workspace = FGDB

    for rstr in arcpy.ListRasters("G*"):
        try:
            newrstr = "G" + "2014" + rstr.replace("g", "")
            arcpy.Rename_management(rstr, newrstr)

        except:
            print "Failed to rename: " + rstr

    arcpy.env.workspace = origdir
    return None
# ----------------------------------------------------------------------------- #

This function works fine for a small to medium amount of rasters within a FGDB and takes about 1.5 seconds per raster renamed. I run into issues when there are a large number of files to be renamed. 1,200 files takes > 10 hours. At the same rate as the smaller raster sets this should take < 1 hour...

Does anyone have an idea on how to make this process go faster?

If not I'll set it up to run over the weekend.

I should mention that I am using ArcGIS 10.2.2.


texts:
I'm trying to use CQL but I guess I do it wrong because I get pink layer in return and error in geoserver :

  Binary geometry filter, but first expression is not a property name?
  (it's a class org.geotools.filter.LiteralExpressionImpl)
  java.lang.IllegalArgumentException: Binary geometry filter, but first
  expression is not a property name?

This is how I try to do the DWITHIN (get all lines at 200 meter from buildings)

wLayer.mergeNewParams({'CQL_FILTER': "DWITHIN('geom',querySingle('Buildings','geom','INCLUDE'),200,meters)"})

What am I doing wrong ?


texts: Using script tool in python add-in wizard?

texts:
I am trying to use a script tool in python add-in wizard to make a toolbar which will have button options. The script tool is working fine in the toolbox and I am using them in arcgis 10.2.2. The script tool I am using doesn't interact with the map. It takes in some input and run some geoprocessing. I want to add the script tool in a toolbar so that I can interact with the output of the geoprocessing. When I added it to the toolbar it does show up in the toolbar but doesn't run the script tool when clicked.

Can anyone suggest how I can use scripted tool in python add-in wizard?


texts: Using Function in Field Calculator gives error 000539?

texts:
I'm trying to populate a field using a function. When I run this in Pyscripter, it works fine:

d = {}

# Function #
def find_dups(val):
    d[val] = d.setdefault(val, -1) + 1
    return d[val]
#

expression = "find_dups(!LINK_ID!)"
arcpy.CalculateField_management(streets, "dups", expression, "PYTHON_9.3", "#")

However, when running it from arc toolbox I get this...

  ExecuteError: ERROR 000539: Error running expression: find_dups(!LINK_ID!)

  NameError: name 'find_dups' is not defined

  Failed to execute (CalculateField).*

So it looks like it can't find the function. I tried using the field calculator in arc to get a python snippet and it gave me this...

arcpy.CalculateField_management(streets,"dups","find_duplicates(!LINK_ID!)","PYTHON_9.3","d = {}/n/ndef find_duplicates(val):/n d[val] = d.setdefault(val, -1) + 1/n    return d[val]/n")

When I tried to put this into my script, however, it didn't work??

  ExecuteError: Failed to execute. Parameters are not valid.

  ERROR 000989: Python syntax error: Parsing error

  SyntaxError: invalid syntax (line 1)


texts: ArcGIS Javascript API: Get Map Object from DIV

texts:
So I am currently writing an application using the ArcGIS JavaScript API, and I want to access the map object that is held inside of my div from a js module I have. I was thinking it would be a simple JQuery call but it doesn't seem to be the case, as when I try to use map.centerAndZoom(), it states that it is not a member of the object that returned by that call.

Is there anyway to be able to reference the map object that was created in my main index.html elsewhere?


texts: Why does this Mapserver layer not overlay correctly onto Google Maps?

texts:
I'm using mapserver to serve an overlay of weather data onto a google map (V2) and the weather data is not displaying where it should.

The image below shows the problem i am having. The black lines are the land border of Alaska, as you can see, they do not line up.  The brown lines represent lat long lines.

As you can see the relationship between the lat lon lines and the borders of Alaska align correctly, however they don't match up with google.

The WMS request i am using looks like the following:

http://mysite.com/cgi-bin/mapserv?map=/opt/fgs/www/htdocs/mymap.map&REQUEST=GetMap&SERVICE=WMS&VERSION=1.1.1&LAYERS=radarALASKA_TEST,Graticule&STYLES=&FORMAT=image/png&BGCOLOR=0xFFFFFF&TRANSPARENT=TRUE&SRS=EPSG:4326&BBOX=0,66.51326044311185,45,79.17133464081945&WIDTH=256&HEIGHT=256&reaspect=true

The bounding box is built using the following javascript:

var lULP = new GPoint(a.x * 256, (a.y + 1) * 256);
var lLRP = new GPoint((a.x + 1) * 256, a.y * 256);

var lUL = G_NORMAL_MAP.getProjection().fromPixelToLatLng(lULP, b, false);
var lLR = G_NORMAL_MAP.getProjection().fromPixelToLatLng(lLRP, b, false);

var lBbox = lUL.x + "," + lUL.y + "," + lLR.x + "," + lLR.y;

and this is my mapfile

MAP
#CONFIG "PROJ_DEBUG" "ON"
#CONFIG "CPL_DEBUG" "ON"
#CONFIG "MS_ERRORFILE" "/opt/fgs/www/htdocs/airlog/logs/ms_error.txt"
#DEBUG 5
IMAGETYPE PNG

SIZE 400 300
SHAPEPATH "/opt/fgs/www/htdocs/data"
#FONTSET   "/opt/fgs/www/htdocs//data/fonts.txt"
IMAGECOLOR 255 255 255
TRANSPARENT OFF

WEB
METADATA
"wms_title"           "WMS Demo Server"
"wms_onlineresource"  "http://10.142.62.153/cgi-bin/mapserv?map=/opt/fgs/www/htdocs/mymap.map&"
"wms_srs"             "epsg:4326 epsg:900913"
END
END

PROJECTION
"proj=merc"
"a=6378137"
"b=6378137"
"lat_ts=0.0"
"lon_0=0.0"
"x_0=0.0"
"y_0=0"
"k=1.0"
"units=m"
"nadgrids=@null"
"wktext"
"no_defs"
END

EXTENT -20037508.34 -20037508.34 20037508.34 20037508.34
UNITS meters

LAYER
TYPE RASTER
STATUS ON
#NAME GISIMAGE_RADAR_CUR_ALASKA_TEST
GROUP radarALASKA_TEST
DATA "/opt/fgs/www/htdocs/Bristow/GISIMAGE_RADAR_CUR_ALASK.png"

PROJECTION
"proj=latlong"
"a=6367451.5"
"b=6367451.5"
END
UNITS DD
PROCESSING "DITHER=YES"

# For Mapserver 3.x, we would use
# OFFSITE 0
OFFSITE 0 0 0
END

SYMBOL
NAME "circle"
TYPE ellipse
FILLED true
POINTS
1 1
END
END

#  Graticule (lines of constant lat-lon)
LAYER
NAME "Graticule"

PROJECTION
"proj=latlong"
"a=6378137"
"b=6378137"
#"ellps=WGS84"
#"datum=WGS84"
END

TYPE LINE
STATUS ON
CLASS
COLOR 180 180 180

LABEL
TYPE BITMAP
SIZE SMALL
COLOR 90 90 90
POSITION AUTO
OUTLINECOLOR 255 255 255
END
END

GRID
#  In lat-lon coordinate system, lines will always be straight (need 2 pts)
MINSUBDIVIDE 2
MAXSUBDIVIDE 180
# MAXINTERVAL must be < MININTERVAL or mapserver will hang.
# MININTERVAL 15
# MAXINTERVAL 10
# MINARCS 18
# MAXARCS 10
LABELFORMAT '%.f'
END
END
END

I think the issue is related to the output WKT in the map file but i'm not really sure and i would certainly appreciate help from those far more experienced than me at this stuff!

Thanks everyone in advance, and if i have neglected to include some critical information, please let me know :)


texts: Transfering GPX files from Basecamp to CSV files for use in ArcMap

texts:
How do I transfer GPX files from Garmin Basecamp to CSV files for use in ArcMap?

I have saved all my points and tracks in Basecamp and now need to plot them on ArcMap.


texts: How to modify combo boxes lists in .esriAddIn, with no access to Visual Studio or SDK

texts:
I have an .esriAddIn project made for ArcGIS 10.0, and I want to modify the combo boxes lists with up-to-date values. My organisation has moved forward with Visual Studio 2013, ArcGIS 10.3, and relevant SDKs. The issue is that the esriAddIn project created by these versions doesn't work in ArcGIS 10.0. ESRI no longer provides older versions (to replicate ArcGIS 10.0 environment). So, here is the question:

How do I modify the combo boxes in the .esriAddIn project so as to be used in ArcGIS 10.0? I have no access to past Visual Studio and relevant SDKs, and the newer ones create esriAddIn that is not accepted by ArcGIS 10.0.


texts: After dissolving lines, old node/vertex order does not reset, causing issues with chainage

texts:
In QGIS, I am trying to dissolve multiple segments of a river feature and then calculate chainage on the resulting single feature.  However, it seems like the internal node/vertex order gets preserved sometimes from the old features, which then causes the chainage plugins (both "QChainage" and "Locate Points Along Line") to calculate the chainage measurements incorrectly along the old segments, rather than understand the new dissolved segment as a single line.  Is there a method to correct or reset the internal vertex ordering (I guess that's what it is... I'm not an under-the-hood GIS expert by any means) after the dissolve so this issue won't occur?

Image below displays issue, this is a single line feature resulting from dissolving 2 line segments. The chainage shows it 'beginning' within the line itself near the upper left corner, then jumping back to the other end.  This issue occurs when dissolving multiple line features representing a single river in the NHD flowline dataset.


texts: Calculating the correlation coefficient of one pixel with the surrounding pixels in same raster

texts: I am working on a problem in which I have 20 years of daily data (raster). I have a point shapefile with 123 point locations. Now first I need to select the particular pixel on which the point location is falling and find the correlation coefficient between that pixel and pixels in the 300km radius area of that pixel. This has to be done for all the 123 point locations on all the raster files of daily data (ignoring no data pixels).


texts: I'm trying to work out on how I can insert or overwrite points that are stored in a features/ Layer? I've created a String/Text document which is already converted into GML/XML ready to be pushed to a WFS. But not sure how I can do this, have been reading into it, can't understand how it works.


texts: PostGIS: how to merge contiguous features sharing same attributes values?

texts:
I have a big PostGIS table storing line features.

I would like to merge features that are both connected to each other and have given attribute values identical. Here is below an example of what I'd like to achieve.

The outmost features (black and blue ones) are not merged because there are not contiguous even if their attribute values are the same
Features green and red are merged because they match both conditions
Feature yellow is left as before.

How would you achieve that with a SQL query?


texts:
I have a path of a raster layer. I want to get name of this layer in GDAL.
In below code I want to get aa.tif

from osgeo import gdal

dataset = gdal.Open("/Data/raster/aa.tif")
print dataset.GetName()

But I received below error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/tmp/tmpq3Peqi.py", line 5, in <module>
    print dataset.GetName()
  File "/usr/lib/python2.7/dist-packages/osgeo/gdal.py", line 1804, in <lambda>
    __getattr__ = lambda self, name: _swig_getattr(self, Dataset, name)
  File "/usr/lib/python2.7/dist-packages/osgeo/gdal.py", line 74, in _swig_getattr
    return _swig_getattr_nondynamic(self, class_type, name, 0)
  File "/usr/lib/python2.7/dist-packages/osgeo/gdal.py", line 69, in _swig_getattr_nondynamic
    return object.__getattr__(self, name)
AttributeError: type object 'object' has no attribute '__getattr__'

How can I do it?


texts:
How can I join these lines depending on the Trigger field. The condition is as follows, the fields must be ordered in ascending order after the fid field from 1 to 68. Thus, I want to join all my lines starting with the first one that has trigger 1 and the next ones that have the value 0 and stop at the first value with 1 and start again.
Edit: I understand. But how can I solve this problem? How to join these lines? Trigger refers to the angle. If it is 0 then the angle between 2 segments is less than x degrees, and if it is 1 the angle is greater than x degrees. So, I want if the first line is not a trigger, I want to join it with the first line that is the trigger.


texts:
I want to create a wrapper for my server. (I use mapserver and apache on windows server)

In this document ( http://mapserver.org/cgi/wrapper.html ) refer to suggestions:

1st :

#!/bin/sh
MAPSERV="/path/to/my/mapserv"
MS_MAPFILE="/path/to/my/mapfile.map" exec ${MAPSERV}

2nd :

#!/bin/sh
MAPSERV="/path/to/my/mapserv"
MAPFILE="/path/to/my/mapfile.map"
if [ "${REQUEST_METHOD}" = "GET" ]; then
  if [ -z "${QUERY_STRING}" ]; then
    QUERY_STRING="map=${MAPFILE}"
  else
    QUERY_STRING="map=${MAPFILE}&${QUERY_STRING}"
  fi
  exec ${MAPSERV}
else
  echo "Sorry, I only understand GET requests."
fi
exit 1
# End of Script

1. Are that run on windows server?

2. in that document say : "Put the script in your web server’s cgi-bin directory, and make it executable". what is the name and extension of saved script ?

also in MapServer WMS Document in url "http://mapserver.org/ogc/wfs_server.html" refer to some suggestions :

Apache ReWrite rules (using Apache mod_rewrite)
Apache environment variables - MS_MAPFILE
Apache SetEnvIf
ASP script (IIS - Microsoft Windows)
Mapscript wrapper
Wrapper script (Unix)

3. Which is best suggestion? why? What are the advantages of each?

4. "1st suggestion" and "Apache mod_rewrite" are seems easy and good Enough. Is not it?


texts: How do I calculate average depth of a 250m square using points at 125m spacing using QGIS

texts:
I have a set of points spaced at 125m intervals across an area, an average depth value is assigned to each point.
I need to calculate from these points a single average depth value for a 250m square, I started doing this manually but there are 220 points, and would take a long time to do.
I was wondering if there is a tool in QGIS that I could use to achieve this?


texts:
How do I bring Ordnance Survey's OS Vector Tile product directly into QGIS via their API (https://osdatahub.os.uk/docs/vts/overview)?
It feels like this should play very easily with the vector tiles layer add, but having tried both the plain URL with authentication separate, or with the key embedded, get errors both with generic and ArcGIS style connections. All the examples they have are for building web apps, but I'm looking at high granularity analysis instead so they don't work for me. (Examples here, tech spec here).


texts: Increasing/allocating RAM for ArcGIS usage?

texts:
In my x64 laptop I have 16G of RAM.
I see arcGIS is using 4G, I have 11G of RAM available.
Also, I was trying to process an 8G raster and I've got a not enough memory error.

Is there anyway I can make more RAM allocate/available for arcGIS usage?


texts: Select features from different vector layers at once?

texts:
I have two vector files in the layers panel in QGIS (V2.8.2, Wien).

I am trying to select points from both vector files simultaneously to view in both attribute tables - both tables are set to show selected features only.

The only way I can do this is to select one layer in the Layers Panel, then select features by rectangle, then select the second layer in the Layers Panel and reselect again using the rectangle.

Is there a way to select features from both layers and only using select by rectangle once?


texts: Adding offset to GeoJSON layer in Leaflet?

texts:
Using the example GeoJSON tutorial code at leafletjs.com.
Change the myLines variable so that the two lines overlap, but run in opposite directions as follows:

var myLines = [{
    "type": "LineString",
    "coordinates": [[-100, 40], [-105, 45], [-110, 55]]
}, {
    "type": "LineString",
    "coordinates": [[-110, 55], [-105, 45], [-100, 40]]
}];

This overlapping is a common situation with lineString shapes, especially roads.  In QGIS/ArcMap, one can apply an Offset = X to the layer such that a positive offset shifts each line feature X (pixels, map units) to the right, in the direction of travel (starting coordinate to ending coordinate), and a negative offset shifts the line to the left, in the direction of travel.

There does not appear to be an option in Leaflet which would allow one to apply a similar offset in the style to these lines so that they are not shown overlapping.  Markers seem to have this option, but not lines.  I was hoping the following code would work, but it doesn't since "offset" is not a valid option.

var myStyle = {
    "color": "#ff7800",
    "weight": 5,
    "opacity": 0.65,
    "offset": 1.5
};

L.geoJson(myLines, {
    style: myStyle
}).addTo(map);

Does this option exist or is it perhaps called something else?


texts: Can I have two distinct workspaces with two different layers reference the same data store in Geoserver?

texts:
Using geoserver 2.7.2 on Ubuntu 14.04 64 bit.

I have a large image pyramid I'd like to use in two workspaces simultaneously. Do I have to create two separate "stacks" of folders, stores, and layers?

or can I create two data stores (one for each workspace) that reference the same folder?


texts: Convert line shapefile to raster, value=total length of lines within cell

texts:
I have a line shapefile representing a road network. I wish to rasterize this data, with the resulting values in the raster showing the total length of the lines that fall within the raster cell.

The data is in the British National Grid projection so the units will be meters.

Ideally I would like to perform this operation using R, and i'm guessing that the rasterize function from the raster package would play a part in achieving this, I just can't work out what the function applied should be.


texts: Splitting shapefile into multiple by value from field in QGIS

texts:
I am working in QGIS. I have a shapefile and its table and I need to export each row as a new file called "Fecha" field:

Finally, I will get 6 shapefiles such as 20130612.shp, 20130714.shp, ...etc


texts: Accessing and/or manipulating table data in QGIS composer label?

texts:
When setting up a map to print using the QGIS Composer, how can I insert a label linked to a specific table record?

I have been using the latest development build (1.9.0) which has the new Expression Builder functionality, but still can't figure out how to do this.

Example: I have a paddock layer with hundreds of paddock areas scattered over the country. Each area polygon is linked to a farm ID attribute. I want to create a composer view for each farm. On each view, I want a label that inserts the total hectares of all paddocks for that one farm.

Either I have to get the Expression Builder to filter then add up areas and produce a total, or do this separately in a database table which is easy enough - but then how to access it?


texts: SLD_BODY in WMS Layer params is not working

texts:
I am trying to set style to my WMS layer in my script like this :

http://opendem.info/ol3/wms_sld.html

but I got nothing.

This is my code :

var sldBody = '<?xml version="1.0" encoding="ISO-8859-1"?><StyledLayerDescriptor version="1.0.0" xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"><NamedLayer><Name>opengeo_arbousiers2</Name> <UserStyle> <Title>A azure polygon style</Title><FeatureTypeStyle> <Rule> <Title>azure polygon</Title><PolygonSymbolizer><Fill><CssParameter name="fill">#0033cc</CssParameter></Fill> </PolygonSymbolizer> </Rule> </FeatureTypeStyle> </UserStyle></NamedLayer></StyledLayerDescriptor>';
var map = new ol.Map({
  layers: [  new ol.layer.Image({
          source: new ol.source.ImageWMS({
            url: 'http://localhost:8080/geoserver/opengeo/wms',
            params: {'LAYERS': 'opengeo:regions1','SLD_BODY' :sldBody},
            serverType: 'geoserver'
          })})],
  target: document.getElementById('map'),
  view: new ol.View({
    center: [-1095791.453557, 3422374.879112],
    maxZoom: 19,
    zoom: 5
  })
});


texts: How to test for correlation between two point data sets

texts:
I want to determine if there is a correlation between weed plant species density and habitation (buildings).
Both data sets are points.  The weed species has an 'abundance value' (i.e., one point may represent any value captured in an attribute field). The habitation point data set is one 'building' per point.
I want to test both; for 1) all species relationship to buildings, and 2) per single species to buildings (i have 115 species). The relationship should be tested up to about 3000 m from each species point.  Locations are not coincident.

I only have ArcGIS 10.4 with Spatial Analyst, Geostatistical Analyst, and other common extensions. I would also be able to use QGis.

I have little knowledge of R, python, etc.

My hypothesis would be that 'weed species are more abundant close to habitation and more so the denser habitation is'.


texts: How to remove colors in an RGB raster to display only contour lines and roads for a basemap

texts:
I have an old topo map which I want to use as a base map and then display transparent polygons of geologic units on top.The green and pink on the base map alters the color of the transparent polygons. I have been able to use the color map on other rasters to select a particular color and make it not display color, however, this basemap is an RGB 3 band raster without a color map when I expand the contents.

Additionally, is there any way for me to isolate the contours and roads with the TIFF I have or do I need to do a conversion to the file before import? If I convert the file, will I lose my referencing data?

I am using ArcMap 10.2.2


texts: Alternatives to Ushahidi for Crowdsourcing?

texts:
We all know that Ushahidi offers a versatile solution for crowdsourcing information and its effective depiction in maps.

I am interested to know if there is any application(free and open source) on par with or superior to Ushahidi!

I would be using the application to collect data from users via tweets, web forms, emails, SMS and map them using the location information provided by them; all of which can be accomplished by Ushahidi.

My deployment would be a permanent online deployment to collect ongoing information and not the one which would be live for few days/weeks.

Hence, one of the major limitation I foresee with Ushahidi is the requirement of manual intervention for approval of reports, which may need more personnel in the event of large number of responses.

Is there any possibility to build a custom solution for such purpose using different open source tools?


texts:
As a part of my PhD I have to map the rivers of Northern Peloponnese, Greece.
I have downloaded two DEMs from Nasa Reverb, which I connected with "mosaic to new raster".
I have started the usual process of "fill", "flow direction" and "flow accumulation" in order to be able to analyze the rivers of the area.
But unfortunately I'm facing two problems.

First one is that the results I get from Flow Direction tool are wrong and I can't edit them, no matter what I try.
As you can see at the screenshot, there are some horizontal flow lines that make no sense.
I tried editing the values (tried to make arcmap exclude values lower than 1) of the original raster and then using "fill" but that doesn't prevent the "fill" tool from ignoring these areas.

Second problem (which seems to be connected to the first one) is that I decided to move to the next step and use "Flow Accumulation" anyway as some areas seemed to be fine and see what happens.

As you can see from the screenshot the flow accumulation lines follow the shoreline shape and then extend in the sea. As far as I can understand, these flow accumulation lines follow the flow direction ones so the problem is created when I use the "flow direction" tool.

Can you please help?
I'm fairly new to this so excuse my ignorance.

If it's a problem with the original data and the problem can't be corrected, is there a way around it?


texts: Generate Service Area tool skipping polygon generation when run with Arcpy command

texts:
I created a tool that uses service areas to group points that are within a certain distance from one another along a network (in my case, a river system). Overall, the tool works. However, when I run it, the Generate Service Areas geoprocessing step does not produce polygons around all input point 'facilities' when run with the Python code.

I tried running the same datasets manually with the Network Analyst toolbar and window and all polygons were successfully created. I also tried running the same dataset through the Generate Service Areas geoprocessing tool and it also generated all polygons as expected.

However, when I run the Python code snippet from the Python window or from within my custom tool, it randomly skips polygon creation for several of the point 'facilities'. This is my code snippet for the generate service areas step:

species_pts = 'species_pts'
sep_distance = 500
scratchGeodatabase = r'C:\Users\mmoore\Documents\ArcGIS\Default.gdb'
arcpy.GenerateServiceAreas_na(species_pts,sep_dist,"Meters",network,os.path.join(scratchGeodatabase,"service_area"),Polygons_for_Multiple_Facilities="MERGE",Detailed_Polygons="DETAILED_POLYS")

Anyone have an idea as to why this may be somewhat randomly skipping polygon creation around several points?


texts:
I have setup MapServer (MS4W) on my local Windows machine, and can access the map with this link http://localhost/cgi-bin/mapserv.exe?template=openlayers&layers=all&map=C:/ms4w/apps/osm/basemaps1/osm-google.map

My requirement is I want directly jump to a lat long with zoom value and show marker on lat long.

I looked into OpenLayers and with a overview I think it is good for my requirement. How can I setup my map=C:/ms4w/apps/osm/basemaps1/osm-google.map with OpenLayers?


texts: Creating buffer around points (layer + map in the same PCS)?

texts:
Currently creating a map of fire stations in LA and want to create a 50-mile buffer around them. My project right now is in NAD 27 / EPSG:26799 (Projected CS). I also managed to save my shapefile with all the points into the same PCS.

However, when I try to create a buffer around these points, I am still getting a huge thing that takes up the whole map. What exactly am I doing wrong here?


texts:
I have a CSV file with a column that is WKB, called WKB.

When I import it to QGIS (using "Add Delimited Text Layer"), there seems to be no option for WKB, and the file imports as pure CSV - no geocoded polygons show up on the map, and it thinks the WKB field is a string.

Is there any way to import WKB to QGIS?


texts: Merge parts of line as one single atlas feature

texts:
While you can build a atlas with one page for every single feature in a coverage layer, it seems to be more difficult, to build one atlas page from a couple of features.

Example: I have 2 tracks and each is split in several parts. In attribute table it looks something like this:

Feature   track
Part 1      1
Part 2      1
Part 3      2
Part 4      2
Part 5      2

If I build a atlas filtered by track, I will get 2 pages for every part of track 1 and 3 pages for every part of track 2.
But I need just one for all parts of track 1 and one for all parts of track 2.
Of course I could merge the parts but I can/will not do any changes to the layer. I could also build a new file with merged parts - but also this won't be a good solution, too. I'm looking for a more "dynamic" way with PostGIS, Python, QGIS


texts:
I have written a Python script for QGIS 2.14.3 which reads in every .txt file in a folder containing hundreds of .txt files (format: "X Y Z" whitespace-separated)and converts them to ESRI shapefiles.

As a next step I want to generate contour lines for every shapefile. Therefore I have installed the Contour plugin (V.1.4.2) from
https://github.com/ccrook/QGIS-Contour-Plugin

How can I call the Plugin from my script?

Can I somehow import my plugin folder to make my installed plugins accessible?

When I look into /homefolder/.qgis2/python/plugins/contour there are five .py files...

Of course I could do the whole workflow by hand with the GUI, but that would take too long.

If somebody needs the script mentioned above (which converts the .txt files to shapefiles), here is the code:

# Defines Input parameters which are displayed in the GUI
##Input_Folder=folder
##Unbuffered_PolygonLayer=vector
##output_Folder=folder

import os
from qgis.core import *
from PyQt4.QtCore import *

indir = Input_Folder
outdir = output_Folder

for dirpath, dirnames, filenames in os.walk(indir):
    for file in filenames:
        if file.endswith('.txt'):
            fullname = os.path.join(dirpath, file).replace('\\', '/')
            filename = os.path.splitext(os.path.basename(fullname))[0]
            #type=(csv|regexp|whitespace|plain) , skipLines=n, useHeader=(yes|no), trimFields=(yes|no), skipEmptyFields=(yes|no), decimalPoint=c, xField=column yField=column [ If the useHeader is no (ie there are no column names), then this is the column number (with the first column as 1).]
            # crs=crsstring, spatialIndex=(yes|no)
            # see: http://qgis.org/api/2.14/classQgsVectorLayer.html

            uri = ( 'file:///%s?type=%s&skipLines=%d&useHeader=%s&trimFields=%s&skipEmptyFields=%s&decimalPoint=%s&xField=%s&yField=%s&crs=%s&spatialIndex=%s' % (fullname,   "whitespace", 0, "Yes", "Yes", "Yes", ".", 1, 2,  "EPSG:31468", "Yes"))
            layer = QgsVectorLayer(uri, "LayerName", "delimitedtext")

        if layer.isValid():
            print "Loaded Layer %s is valid..."%(fullname)

            #QgsMapLayerRegistry.instance().addMapLayer(layer)
            #QgsVectorFileWriter.writeAsVectorFormat(layer, outdir + '/' + filename + '.shp', 'CP1250', None, 'ESRI Shapefile')

        else:
            print "Failed to load Layer %s" % (fullname)


texts:
Is there an archive through which I can access older versions of SoilGrids? I'm specifically looking for one representative dataset for each year since 2015.
Also, is SoilGrids updated on a regular basis, and if so, how often is it updated?


texts: Trying to place a Point geometry on OSM map, unsuccessfully

texts:
I have a piece of code that takes a map layer and places a vector layer on top of that.
Currently, I'm just trying to place a few points on top of the map, unsuccessfully.

I managed to place 1 point but it would only stay in the middle of the map.

I took all the html code and simplified it.
Here is the code for the whole page, it's standalone, so you can just copy paste into a new html page and see the problem.

To reiterate, how do you place the point at some location on the map? (I have lon/lat)

Thanks!

Update

Just to make it easier for readers in the future.
The answer which is marked as the solution below only hints at the full solution.
The 'full solution' is to modify the point creation code as such:

// Make Point
coordinates.push(new OpenLayers.Geometry.Point(33, 33).transform(new OpenLayers.Projection("EPSG:4326"), map
                        .getProjectionObject()));


texts: What is a "feature" in Getis-Ord G* hot spot analysis?

texts:
From the official site, it says "The Hot Spot Analysis tool calculates the Getis-Ord Gi* statistic (pronounced G-i-star) for each feature in a dataset" and "xj is the attribute value for feature j".

Say the entire map is divided into 100 smaller grid units, but only 50 of these 100 units have an input value. Do I have 100 or 50 "features"? For the units that don't have an input value, do they also go through the G* statistics? If so, does it mean it automatically assumes the value of 0?


texts: Performing SDNA 3d network analysis with ArcScene

texts:
I am trying to use SDNA commands to run the 3d network analysis on the ArcScene Shape file. I can view the shape file as 3d, however when I run the SDNA Prepare command, the below error message comes.
Is there anything I am missing or doing wrong?
The same process works perfectly fine for 2d files in ArcMap.


texts: Creating time series (NDVI) for each polygon contained in feature collection in Google Earth Engine?

texts:
I have an image collection of Landsat 5 and a feature collection of thousands polygons and points (in the example, I only have three). I would like to compute the NDVI (and NDWI) mean for each polygon over time like if I was using: ui.Chart.image.series;

If I am testing only one polygon with this code:

var p1= ee.Geometry.Point([-78.55995626672507,35.05443673532838])
var pol = ee.Geometry.Polygon([[[-78.57239414626946,35.01247143741747],
[-78.57186843330254,35.012559309453266],
[-78.57199717933526,35.01277020195395],
[-78.57253362113823,35.01272626606113],
[-78.57239414626946,35.01247143741747]
]]);

var ens = [
ee.Feature(pol, {name: 'Thiessen'})
];

var col =  ee.FeatureCollection(ens)
print(col)

// NDVI: B4 and B3
var addNDVI = function(image) {
var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI');
return image.addBands(ndvi);
};

// Apply the cloud mask and NDVI function to Landsat 5 imagery and   print the chart
 var l5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_TOA")
      .filter(ee.Filter.calendarRange(1985,2007,'year'))
      .filter(ee.Filter.calendarRange(1,1,'month'))
      .filterBounds(p1)
      .map(addNDVI)

print(ui.Chart.image.series(l5.select('NDVI'), col, ee.Reducer.mean(), 30));

With this code I obtained this figure. I would like to obtain the same type of figure with several polygon (one figure contains all the time series of the polygons).

I tried this code:

var p1= ee.Geometry.Point([-78.55995626672507,35.05443673532838])
var p2= ee.Geometry.Point([-78.5725093420931,35.05908805245044])
var pol = ee.Geometry.Polygon([[[-78.57239414626946,35.01247143741747],
[-78.57186843330254,35.012559309453266],
[-78.57199717933526,35.01277020195395],
[-78.57253362113823,35.01272626606113],
[-78.57239414626946,35.01247143741747]
]]);

var ens = [
ee.Feature(p2, {name: 'Thiessen'}),
ee.Feature(pol, {name: 'Thiessen'})
];

var col =  ee.FeatureCollection(ens)
print(col)

// NDVI: B4 and B3
var addNDVI = function(image) {
var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI');
return image.addBands(ndvi);
};

// Apply the cloud mask and NDVI function to Landsat 5 imagery and    print the chart
var l5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_TOA")
      .filter(ee.Filter.calendarRange(1985,2007,'year'))
      .filter(ee.Filter.calendarRange(1,1,'month'))
      .filterBounds(p1)
      .map(addNDVI)

//Create a graph of the time-series.
var graph = ui.Chart.image.seriesByRegion({
imageCollection: l5,
regions: col,
reducer: ee.Reducer.mean(),
scale: 30,
})
print(graph)

This code provides this figure

The last code presents the graph as I would like. However, it does not compute what I want. For the polygon pol, I should have the same graphs with the two codes which is not the case. How could I have the same computation made with the code 1 but presented as the code 2?


texts: With a free CartoDB account, visitors can leave comments on your map/data.  Is there any way to turn this feature off?


texts: Using Select Layer By Attribute in ArcGIS ModelBuilder multiple times?

texts:
I want to use Select Layer By Attribute in ArcGIS ModelBuilder multiple times.
It works everything fine until a new feature is created (e.g. after Dissolve or Multipart to Singlepart). Then I can't apply Select Layer By Attribute anymore. The Layer Name or Table View of the Dissolved feature is not selectabled in the Select Layer By Attribute menu. My workaround at the moment is to build many little models.

Question: How can I use Select Layer By Attribute in ModelBuilder after a new feature is created?


texts: QueryTask doesn't return anything ArcGIs Javascript 3.24

texts:
I have a query that I believe I am setting up correctly, but maybe I'm missing a vital property. Does the where clause need to be set? I want to query all features, not a subset of them, so I don't really need a where restriction.

Every time I run this query it goes to the error callback where the message given is unhelpful.

printQLyr = new QueryTask("https://xxx/MapServer/3");
printQuery = new Query();
printQuery.returnGeometry = false;
printQuery.outFields = ["Name"];
printQLyr.on("complete", printQResults);
printQLyr.on("error", errorResults);
printQLyr.execute(printQuery);

function printQResults(featureSet) {
    var resultFeatures = featureSet.features;
    console.log(resultFeatures[0]);
}

function errorResults(message) {
    console.log(message);
}

What I get back in the console is this:

Error at Object.A.m.load (https://js.arcgis.com/3.24/init.js:1028:392) at https://js.arcgis.com/3.24/init.js:87:429 at c (https://js.arcgis.com/3.24/init.js:103:393) at d (https://js.arcgis.com/3.24/init.js:103:182) at b.Deferred.resolve.callback (https://js.arcgis.com/3.24/init.js:105:10) at c (https://js.arcgis.com/3.24/init.js:104:96) at d (https://js.arcgis.com/3.24/init.js:103:182) at b.Deferred.resolve.callback (https://js.arcgis.com/3.24/init.js:105:10) at https://js.arcgis.com/3.24/init.js:92:118 at h (https://js.arcgis.com/3.24/init.js:107:277)


texts: Transformation of geometry on select Postgresql

texts:
I am trying to run a transformation of a geom via a select query, the layer has been digitised in WGS 4326 and when I run my select to filter the table by my criteria, I would like to transform the CRS from WGS84 (4326) to UTM36S (32736)

My final Query is designed make a new table of attributes from the first changing the CRS 4326 - 32736

I have two steps I cannot seem to join

SELECT ST_AsText(geom) FROM foo  which provides the list of geoms as text
MULTIPOLYGON(((32.841274 -25.151649,.............32.841274 -25.151649)))

And with the resultant Text entered in:

SELECT ST_AsText(ST_Transform(ST_GeomFromText('MULTIPOLYGON(((32.841274 -25.151649,.............32.841274 -25.151649)))',4326),32736)) As wgs_geom;

I get the result:

MULTIPOLYGON(((484002.707664207 ............... 7218250.78539813)))

But:

SELECT ST_AsText(ST_Transform(ST_GeomFromText FROM (SELECT ST_AsText(geom) FROM foo) ,4326),32736)) As wgs_geom;

or

SELECT ST_AsText FROM a (ST_Transform(ST_GeomFromText FROM (SELECT ST_AsText(geom) FROM foo) as a ,4326),32736)) As wgs_geom;

doesn't run I've spent a bit of time looking at the syntax of the SQL but can't get it to run...

Any Ideas?

Note: the MULTIPOLYGON are quite complex hence ..........


texts: Seeking Open source Tools for GML in ArcGIS and gvSIG?

texts:
I’m searching open source tools for ArcGIS Desktop 10.5 and gvSIG 2.4 to load and write complex GML data for example INSPIRE GML.

I have not the possibility to install interoperability or INSPIRE extension for ArcGIS and I’m looking for other solutions.

gvSIG instead read GML file but doesn’t write them, I found the extension for INSPIRE cadaster data but I need something for more general topics.


texts: Cannot convert string to int using ArcPy Calculate Field?

texts:
When I load a csv file to join with ArcPy, it converts all string fields that are numeric to integers. Therefore, I have to add another field within the feature class that I'm using for joining that is an integer version of the one I originally intended to use.

There are two approaches I've tried to convert the contents of the original field to an int:

Approach 1:

arcpy.AddField_management(layerName, "BG10_INT", "LONG", 12, "", "", "", "NULLABLE", "")
expression = "toInt(!BG10!)"
codeblock = """def toInt(value):
    if value.isdigit():
        return int(value)
    else:
        return 0"""
arcpy.CalculateField_management (layerName, "BG10_INT", expression, "PYTHON", codeblock)

Approach 2:

for row in rows:
    if row.BG10.isdigit():
        row.BG10_INT = int(row.BG10)
    else:
        row.BG10_INT = 0
    rows.updateRow(row)

With either approach, I get The value type is incompatible with the field type. [BG10_INT]. How could something so simple as converting string to int be such a challenge?

After several comments and an answer suggesting I define the schema.ini file for the csv file instead of the above conversions, I tried the following approach:

Approach 3

env.workspace = "C:/data/PCA"
f = open("C:/data/PCA/schema.ini", "wb")
text = ["[bg_to_tr.csv]", "ColNameHeader=True", "Format=CSVDelimited", "Col1=BG10 Text Width 12", "Col2=PrefBUC1 Text Width 10"]
for row in text:
    f.writerow(row)

f.close()

env.qualifiedFieldNames = False
inFeatures = "C:/geodatabases/Intermediate.gdb/BG_sample"
layerName = "BG_sample_lyr"
joinTable = "bg_to_tr.csv"
joinField = "BG10"
arcpy.MakeFeatureLayer_management (inFeatures,  layerName)
outFeature = "C:/geodatabases/Intermediate.gdb/BG_sample_w_tr"
arcpy.AddJoin_management(layerName, joinField, joinTable, joinField)
arcpy.CopyFeatures_management(layerName, outFeature)

This approach only works when I join manually (sans ArcPy) within ArcMap. However, for my purpose, this needs to be done with ArcPy. Is there an explanation why this only seems to work with ArcMap? The csv file was joined without adding it directly into the ArcMap TOC.

Here is a sample of the contents of the csv file:

"BG10","PrefBUC1"
"170312106022","0002804408"
"170312106023","0002804408"
"170312105012","0002804408"

If I open the csv file in ArcMap, it recognizes schema.ini and indicates the BG10 column as a String of length 12. When I look at the join output in BG_sample_w_tr, however, it indicates that it read it as data type Double.

If it matters, neither table contains null values for the joined fields.


texts: Showing OpenStreetMap intersections in ArcMap?

texts:
I'm trying to carry out a few tests on OpenStreetMap data to test its quality, using ArcGIS.

One of the things I'd like to investigate involves checking that all road intersections are represented by a node, rather than two lines crossing over each other, with no node.

I've downloaded OSM data from GeoFabrik, but the intersections don't seem to be marked.

Does anyone have any ideas?


texts: How to extract a rectangle corner coordinate in Python script

texts: I want to extract a rectangle corner coordinates of a shapefile in Python script... please suggest me a solution.


texts: Improving performance of shapely crosses

texts:
I am unsure if this is the fastest way to test if LineStrings cross each other. Is there a better method.

for line_1, line_2, in itertools.combinations(lines, 2):
    crosses = line_1.crosses(line_2)
    if crosses is True:
         print(line_1, line_2)

Obviously this gets pretty slow with large files, I have found maybe loading it into an rtree index may speed this up.


texts:
I have Contour lines shapefile which I imported that into GRASS and converted to a DEM using the IDW interpolation technique (r.surf.idw with npoints 250). Now I want watershed analysis to be done where I can get Flow accumulation, Flow direction, basin, Stream order, drainage pattern etc.

When I Try to run r.fill.dir after a while GRASS crashes and the process doesn't give any output. I have tried other hydrologic modelling under raster but it does the same.

Can anybody help me out with the procedure of doing this analysis in GRASS GIS.


texts: Running a Focal Mean or Low Pass Filter or 3 by 3 Mean Kernel in GDAL \ gdal_calc

texts:
I am looking to run the equivalent of Focal Statistics in gdal_calc.  Is it possible?

It is a simple focal I require. I need to run a 3 by 3 kernel over the data (mean) and if possible only on NoData (but if on all that is okay as well).

I can do this in other software but specifically asking if a GDAL_calc solution exists.

This is the exact thing I am trying to do on a singular raster

Con(IsNull("C:\temp\SST_Final_C_NoData.tiff"), FocalStatistics("C:\temp\SST_Final_C_NoData.tiff", NbrRectangle(5,5, "CELL"), "MEAN"), "C:\temp\SST_Final_C_NoData.tiff")

GDAL 1.11.0 and Python 3.4.1


texts: I asked a question a while ago to figure out how to split a line layer with a points layer using QGIS, but cannot, however, get the accepted answer to work due to limitations of my knowledge and the fact that i'm working on Windows 7. I noticed this feature request that dates from about a year ago, requesting a tool that accomplishes exactly what I need. I've been looking for it, don't see it, and am wondering if it is currently an available tool and i'm just not looking in the right spot or not.


texts: Connect multipart polygon with lines to one label

texts:
I have several one-part and multipart polygons that I want to label.

I Set X and Y values to move the label to a certain position and want to connect the label to the polygon using the following expression in a geometry generator:

  make_line(   point_on_surface(  $geometry ),  make_point( "Label_X" ,
  "Label_Y" ))

I want the line to connect from the label to every part of a multipart Polygon.

Is this possible?


texts: Exporting SHP files into AutoCAD readable files?

texts:
I have a full water supply network drawn in a shapefile which I created using Qgis 2.2.0.
Now I need to export it into a format that enables me to open it in AutoCAD.
I am trying to export it as a dxf file, but when I upload it in AutoCAD it shows an empty drawing (although the dxf file is 1.75MB).
The shapefile is a line shapefile with several lines with a different main attribute (the diameter).

Maybe I have to save different dxf for each diameter?


texts: Python syntax error from raster calculator in ArcGIS ModelBuilder?

texts:
I am experiencing a problem with the raster calculator in the ModelBuilder under ArcMap 10.

I need to combine two layers:

a drainage map (issued from a series of previous operations on a DEM) and
a slope map reclassed in 8 classes.

The formula used in the calculator is the following:

Con(("surf_drainee" > 12) & ("surf_drainee" <= 80) & ("pente_rclass" <= 24),"pente_rclass"+1, Con(("surf_drainee" > 20) & ("pente_rclass" < 24),"pente_rclass" + 2, Con(("surf_drainee" > 80) & ("pente_rclass" == 24), "pente_rclass" + 1, "pente_rclass")))

When I enter this in the calculator, using the buttons to enter the expression, it sends back a python syntax error. I discussed it with a colleague, who told me to write it out completely with the keyboard.

This works in the Raster Calculator, when used alone from the Toolbox. But when this is done in the ModelBuilder, it does not connect my two layers to the calculator, but creates numerous duplicates in entry to the calculator (pente_rclass (2), pente_rclass (3)...).

Is there a way of getting around this problem? I was thinking maybe a script could replace the raster calculator, but I don't really have the knowledge to set that up without a few hints.


texts:
Background

Looking to leverage an existing GIS system for exploring organizational data.

Architecture

The following figure represents a high-level overview of the system's desired features:

The most basic usage would be as follows:

The user visits a web site.
The system presents a map (having regions, cities, and buildings).
The user drills-down on the map to a particular building.
The system provides a basic CRUD interface.
The user can view and modify information about personnel (e.g., their assigned teams), equipment (e.g., network appliances), applications, and the building itself (e.g., contact and phone numbers).

Ideally, all the components should be open-source (or otherwise free).

Problem

This must be a small project that needs a quick (but functional) prototype, mostly to confirm whether or not such a system would be useful in the long term.

Questions

What software components would you use to quickly develop a working prototype?
What open-source solutions already exist, if any?

Ideas

Here is what I am thinking:

PostGIS - Define the regions, cities, and sites
Google Maps - Display an interactive, clickable map
geoJSON - Protocol between PostGIS and Google Maps
Seam - CRUD interface

Custom Development

For example, this would entail:

Installation and configuration

Configure SSH for remote logins
Subversion (or git)
PostgreSQL
PostGIS
Java
Tomcat
Seam
JasperReports

Enter GIS information into PostGIS
Aggregate data sources into PostgreSQL database
Develop starting page for map interface
Develop clickable Google Maps interface
Develop summary reports
Develop CRUD interface using Seam for data maintenance

Surely something like this already exists?

Thank you!


texts: Obtaining envelope of data frame in Layout using ArcObjects?

texts:
This is to expand on Transforming layout page units to map units using ArcObjects?

How can i programatically obtain the corner coordinates (in page units) or the envelope of the data frame in layout view?

I don't seem to be able to do this with iActiveView:Extent as it is returning the extent of the data view in page units. The below screen capture shows what I would like to get. The coordinates are roughly; upper left 0.9,10 and lower right 8, 0.9. Even in page layout iActive:extent returns: upper left -3.4, 12 and lower right 14,-3 which is the extent of the entire application window not the data frame as it appears in layout.

Here is what i need:


texts: Switching leaflet image overlay layers shows white blinks

texts:
I am switching between image overlay layers upon user entry (through a slider).

Every time users get into this loop, a fast blink (with white background) is displaying, meaning that the old layer is removed before the new one finished rendering. I was relying on the "load" event, but it seems not affordable. It looks like the "load" event (fired by the added overlay image) is fired before the load completes.

mymap.addLayer(overlay[newlyr]);
overlay[newlyr].once("load",function() {
    mymap.removeLayer(overlay[oldlyr]);
    });

In other applications I have verified that the "load" event is very reliable for tile Layers, but it isn't at all for Image Overlays.
I hope I am doing something wrong.

Any tip?


texts: QGIS raster style rendering incorrectly/unexpectedly

texts:
I have a raster of flow accumulation. To style it, when I set manual values for raster breaks in qgis 2.14 it isn't rendering all values below the value I set in the same colour. My understanding is that all cells below a value level take that colour; however, here the values of the cell remain uncoloured even though they fall below the value I've put in.

When I check the values of the pixels that aren't rendering in colour they are below the value I used. For example, I set the value to 5.13e+11, while the cell in the figure below has a value of 8.19e+9. Shouldn't the values in the high accumulation areas also then become pink (yes, apologies for the hideous pink colour...).

Any ideas where I'm going wrong or what the problem is?

I've tried the same styling in QGIS 2.10 and it works as expected. I assume it's a bug in QGIS 2.14 and have submitted a report (http://hub.qgis.org/issues/14689).


texts: Adding Texture to the sides of an extruded polygon in ArcScene

texts:
I've been using ArcScene 10.2.2 to create simple extrusions of building footprints and now would like to add textures to them to make them more realistic.

I want to add textures to the sides of the extruded polygons but it only adds texture to the top of the polygon and then the sides become stretched versions of that texture.

I've tried and failed many methods of doing this but the only way I found it to convert the polygons to lines and then symbolize the lines with a "3D texture line" then "extruding" using the line width and changing the vertical orientation. This is a very long workaround way to do it.

I don't have access to any other software options (sketchup etc) so was hoping there was a way in Arcscene to do this.


texts: How do I interpret this proj4string and why doesn't changing units change its resuls?

texts:
I'm new to GIS - I have a background in statistics. I have the following proj4string:

'+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0'

Can someone explain:

What the parameters mean (I think lcc means 'lambert conic conformal')
Where I should go to learn how this works/what projections are available?
The input lat, lon coordinates are in decimal degrees. What are the output units? The +units=us-ft suggests feet
I want to focus my analyses on Manhattan, Brooklyn, Queens, the Bronx. As part of my analysis, I want to split my projected map into a grid. Ideally, each cell will be, say, 1/10th a mile square. What projection should I use for this sort of work? Is LCC acceptable?


texts:
I am adding a feature layer of ol.geom.Point geometry as follows:

var location = result.geometry.location;
var lat = location.lat.apply();
var lng = location.lng.apply();
var locationCoordinates = [lat, lng];
var vectorSource = new ol.source.Vector({});
var iconFeature = new ol.Feature({ geometry: new ol.geom.Point(ol.proj.transform(locationCoordinates, EPSG_3857, EPSG_4326)) });
vectorSource.addFeature(iconFeature);
var vectorLayer = new ol.layer.Vector({ source: vectorSource });
map.addLayer(vectorLayer);
vectorLayer.setVisible(true);
var extent = vectorLayer.getSource().getExtent();
map.getView().fitExtent(extent, map.getSize());

But the layer is misplaced, not exactly on longitude and latitude.

I am using EPSG:4326 projection, transforming Google Maps that's supposed to use EPSG:3857 projection.


texts: Try/Except Statement gives different Results in QGIS Console vs. QGIS Tool

texts:
I want my script to do different things depending on whether a particular layer DOES or DOES NOT exist in the Map Layers list.  I thought the try: / except: statements would be the best way to achieve this.

Here is a simplified version of my code:

 checker = ' '

 try:
     vl = QgsMapLayerRegistry.instance().mapLayersByName('mapLayerX')[0]
     boolean = qgis.utils.iface.legendInterface().setCurrentLayer(vl)   # if 'vl' exists boolean == "True"
     checker+='A, '
 except:
     boolean = False
     checker+='B, '

 if boolean == True:
     checker+='C, '

 print checker

I'll look separately at the 2 possible scenarios (mapLayerX absent vs. mapLayerX present).

LAYER IS ABSENT:

Both the Console and the Tool print the same CORRECT result for 'checker': "B, " indicating ONLY the EXCEPT statement is executed (as it should be in those cases).

LAYER IS PRESENT:

FROM THE CONSOLE the code again works perfectly.  The TRY statement executes as it should, setting boolean to 'TRUE' in response to the fact the setCurrentLayer function found "mapLayerX".  This means the IF statement also executes, and the EXCEPT statement is ignored.  The end result: 'checker' is printed as "A, C, "

But run AS A TOOL the code does the OPPOSITE!  The TRY statement appears to be completely ignored somehow (which of course means the IF statement is too) and instead the EXCEPT statement is executed.  The end result: 'checker' is printed as "B, "

 SUMMARY TABLE (results of "print checker"):

                  Console:       Tool:
 Layer Present:   A,C,           B,
 Layer Absent:    B,             B,

*** when the Layer is Present the results should ALWAYS be "A, C, "

Any ideas why the same script runs incorrectly as a Tool (perhaps there are more appropriate Python statements to use here)?  I'm using Q version 2.8.2.


texts: Catching arcgisscripting.ExecuteError and error details?

texts:
I am with ArcGIS 10.1. In my Python script, I need to catch arcgisscripting.ExecuteError and other general Python errors. When a arcgisscripting.ExecuteError is thrown, I also need to retrieve error details.

For example

try:
      # UploadServiceDefinition_server would throw a arcgisscripting.ExecuteError error if an error occurs
      arcpy.UploadServiceDefinition_server(sd, 'My Hosted Services', service,
                                                 '', 'FROM_SERVICE_DEFINITION', 'EXISTING',
                                                 'STARTED', 'OVERRIDE_DEFINITION', 'SHARE_ONLINE',
                                                 'PUBLIC', 'SHARE_ORGANIZATION', '')

except arcgisscripting.ExecuteError:

except:
  general errors
finally:
  do other stuff


texts:
I have the data of the volume per area for Madrid.
from mpl_toolkits.axes_grid1 import make_axes_locatable
def setColorbar(ax, colormap, vmin, vmax):
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    sm = plt.cm.ScalarMappable(cmap=colormap, norm=plt.Normalize(vmin=vmin, vmax=vmax))
    # fake up the array of the scalar mappable. Urgh...
    sm._A = []
    return plt.colorbar(sm, cax=cax)

vmin=0
vmax=70
cmap='jet'
f,ax=plt.subplots(figsize=(10,10))
dfMd.plot(column='volume', alpha=0.5,ax=ax, cmap=cmap, vmin=vmin, vmax=vmax)
setColorbar(ax, cmap, vmin=vmin, vmax=vmax)

I would like to find hotspot of the building density and assign to each polygon a sort of kernel density class associated to the hotspot like this example
What I have done is to define points of interest based on percentiles in this way.
pcts95 = np.percentile(dfMd['volume'], 95)
dfMd95 = dfMd[dfMd['volume'] >= pcts95]
dfMd95 = dfMd95.drop_duplicates(["geometry"]).reset_index(drop=True)
dfMd95 = gpd.geoseries.GeoSeries([geom for geom in dfMd95.unary_union.geoms])
dfMd95=gpd.GeoDataFrame(dfMd95)
dfMd95.columns = ['geometry']
X=[]
Y=[]
for i in dfMd95.index:
  X.append(dfMd95['geometry'].centroid.x)
  Y.append(dfMd95['geometry'].centroid.y)

vmin=0
vmax=70
cmap='jet'
f,ax=plt.subplots(figsize=(10,10))
dfMd.plot(column='volume', alpha=0.1,ax=ax, cmap=cmap, vmin=vmin, vmax=vmax)
setColorbar(ax, cmap, vmin=vmin, vmax=vmax)
ax.scatter(X,Y, color='black', s=75)

However I would like to know if there is a way to estimate kernel density based polygons value as this example


texts: Limiting the number of layers a user can select - Leaflet and WMS

texts: I have a map with around 20 WMS layers but I want the person exploring the map to only be able to have 3 layers active at a time.  For my layer control I am using Styled Layer Control.  What is the best way to do this?


texts: Difference between Identify task & Attribute Info in ArcGIS API for JavaScript?

texts:
I am confused about the difference between Identify task and attribute info task.

According to research, both are displaying attribute information

Interestingly if one point contains multiple attribute then Identify task shows multiple records while Attribute window display only single record
Identify task displays the information if record is hidden in certain cases while Attribute window not.
If you are using Dynamic layers in ArcGIS Java Script API (ver. 2.x & 3.x) then you cannot use Attribute window option while identify task works for both (i.e. feature layer and Dynamic layer)

So what exactly difference between them?  & which is recommended?


texts:
The attached script was created for the geojson file I created that works properly.

<script>
     $(document).ready(function() {
      $(pz).click(function() {
        L.mapbox.featureLayer('map.geojson').addTo(map)
    })
    });
</script>

My question is how can I hide the displayed layer with the click feature?

I would like to ask you some help to create the script.

(I use mapbox.js not mapbox-gl.js)


texts:
I'd like a Linux program which displays ASCII ArcGrid files (NxM rectangular ASCII data grids which may have either integer or floating point values). There is no practical upper or lower limit to these values.

I would like suggestions of a program which could read in this rectangular grid and display it in colour. The program must also have the capability to rescale the data so that I can isolate ranges of interest, it would be best if the rescaling could be nonlinear. (A nice, intuitive interface for this would be bestest.)

This is somewhat of a GIS application, so answers dependent upon conversion to the GeoTIFF format are acceptable.

Gimp did not work for me, so if you suggest it, you should be very specific about how to get it to display the GeoTIFF format (or a rectangular grid)

Thanks!


texts:
I have a shapefile composed of polygons (sample here), many of which have holes in them. I cleaned the polygon layer in GRASS. After cleaning, QGIS's 'Check Geometry Validity' tool returned no errors. I then loaded the shapefile to R as follows:

polys <- readOGR(dsn=path.expand("U:/FandR"), layer="polygon")
polysCRS<- spTransform(polygons, CRS("+init=epsg:32717"))
polysCRS@data$id=rownames(polysCRS@data)
polysCRS.points = fortify(polysCRS, region="id")
polysCRS.df = join(polynsCRS.points, polysCRS@data, by="id")

Unfortunately, upon plotting, I can see that the holes in the polygons are filled.

#plot the polygons!
ggplot(polysCRS.df)+
  aes(long,lat,group=group)+
  geom_polygon(fill="#006600")+
  geom_path(color="white") +
  coord_equal()

Here is an image of a portion of my polygon plot, zoomed in to see the worst portion of the plot.

How can I correct this mistake in my loaded SpatialPolygonsDataFrame?


texts: Is it possible to import an autocad feature to a feature layer without exporting to .shp/.dwg

texts:
Does anyone know if it is possible to import a feature (point, line or polygon) drawn in Autocad directly to a feature layer in ArcMap, so in essence a "feature sucker", rather then exporting the feature to a shp or dwg then importing into ArcMap then copying the feature into the correct feature class.

The datatype that I want to target  is a feature class on ArcSDE (9.2).

Versions are ArcGIS 9.3.1 and AutoCAD 2011.


texts:
I have a CSV edge list such as follows:

Lat1,Lon1,Lat2,Lon2
93.41823,-104.324324,93.41828,-104.324320
93.40110,-104.310552,93.40108,-104.310557
...

Is there a way to import such a CSV file into QGis such that the edges show up as lines?


texts: ESRI JS API - Getting a Screenshot with pure Javascript

texts:
I have a requirement to be able to produce a pdf of the application's map (Esri JS API 3.7).  The map will contain ArcGISDynamicMapServiceLayers, user drawn graphics including text symbols, and possibly a pop-up.

The PrintTask made available at ArcGIS Server 10.1 gets be about 80% of the way there.  However, it does not handle text symbols or popups as far as I can tell.  If it were just the text symbols I had to worry about, I think I would just hack together something with the arcpy mapping module to do the job.  But the popup is much more challenging

I've explored things like BrowserShots, wkhtmltoimage...etc and like with no luck.  I'm wondering if anyone has been able to accomplish something similar.


texts: Stop grid coordinates showing outside map edge

texts:
I'm sure someone else has been looking for a solution to this but can't locate anything. So here goes..
How do I stop my grid coordinates (in the red box) showing outside the map canvas area?


texts:
I have a file of points with Easting and Northing. I have a 1km grid table. Both have the same projection and I can see in a map where the points are layered over the grid squares.

What I need is to find the Grid Square ID that the points fall into.

I have tried all sort of things, nothing I have found in forums/google or QGIS seems to be what I need.

I'm using QGIS Las Palmas.


texts: QGIS graphical modeler add field: error executing algorithm?

texts:
I have to add more than 20 fields to a lot of shapefiles using QGIS 2.4. I used to do it with the tablemanager plugin.  Beside it is time consuming, misstakes can be made. So I want to automate that procedure through the graphical modeler.

I tested it with adding 2 fields and everything worked fine.  But when I added all the algorithms it gave an error:

  Error executing algorithm 2 [Errno 9] Bad file descriptor See log for more details

The log says:

Uncaught error while executing algorithm
Traceback (most recent call last):
Traceback (most recent call last):
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py", line 217, in execute
    self.processAlgorithm(progress)
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\qgis\AddTableField.py", line 75, in processAlgorithm
    self.getParameterValue(self.INPUT_LAYER))
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\tools\dataobjects.py", line 215, in getObjectFromUri
    print _loadedLayers
IOError: [Errno 9] Bad file descriptor

After a few hours I tried to run the model again without changing anything and then this error appeared: Error executing algorithm 6 [Errno 9] Bad file descriptor See log for more details

In the log the same explanation as above.

Can anyone help me understand what's happening here?


texts: Troubles with addFeature in Openlayers by using SelectFeature control

texts:
I'm having troubles with a highlight layer which I created on a map dynamically to put the selected feature on the top of all layers, so I can see it despite of a polygon layer which has been set visible on the SelectFeature method. It is working well so long I don't try to unselect it.
The functions are like:

 var HighlightLayer = new OpenLayers.Layer.Vector('HighlightLayer', {
        visibility: false
        }
    );

function selectedFeaturesSchulenEinzug(event) {
        HighlightLayer.addFeatures([event.feature]);

        // a try with addFeatures but without event.object works well
        //HighlightLayer.addFeatures([OrigVectorLayer.features[0]]);

        HighlightLayer.setVisibility(true);

        //switch on the polygon layer
        PolygLayer.setVisibility(true);
}

  function unselectFeatureEinzug() {

        // switch off the polygon layer
       PolygLayer.setVisibility(false);

       // remove the features and switch off the highlight layer
      HighlightLayer.removeAllFeatures();
      HighlightLayer.setVisibility(false);
  }

It looks like the event.feature object makes the trouble when I click out the layer. Maybe it tries to put a new event.feature into the HighlightLayer and simultaneously removeAllFeatures. However, when I clear the unselect function and it contains nothing, I have the same problem: it tries to do something and after a while the message prompts that the script doesn't answer.
When I try to make the same with a not dynamic feature (like in the script above in the comment), I don't have any troubles and the unselect method can be executed.

Does anybody has experience with a trouble like that?


texts: Why are topocentric coordinates of point A relative to point B not strictly the negative of the topocentric coordinates of point B relative to point A

texts:
I just tried to play around with some geodesy lately, using CartConvert in bash:
$ echo -43.4242920E -22.9145118N 0 | CartConvert -l -43.7954220E -22.7454621N 0
38072.841396 -18768.489365 -141.336023

But if I if I exchange the coordinates given in CartConvert, the result is not exactly the opposite:
$ echo -43.7954220E -22.7454621N 0 | CartConvert -l -43.4242920E -22.9145118N 0
-38119.890905 18672.743845 -141.335023

If we reverse the sign of the two first topocentric coordinates and make a subtraction, some differences appear (here in Python):
import numpy as np
print(np.array([38119.890905, -18672.743845, -141.335023]) -
      np.array([38072.841396, -18768.489365, -141.336023]))
[4.7049509e+01 9.5745520e+01 1.0000000e-03]

47 and 96 meters eastwards and northwards respectively.
Why?


texts: Percentage of area covered by ocean within certain distance using ArcGIS Desktop?

texts:
Here is what I am trying to do using ArcGIS Desktop:

I am studying the coastal area, with several headlands and bays. The goal is to find the best suitable areas to install bee-yards (apiaries). One of the criterion I want to use is the amount of usable forage area that is lost by ocean cover, considering a flight range of 2000 m.

Here some abstract examples to help you visualize what I am trying to explain:

If I would install a bee-yard directly on straight coast I would have a loss of 50% percent of usable forage area.

If I would install a bee-yard in the corner of an bay with an 90º angle I would have a loss of 25% percent of usable forage area.

A 90º headland would be 75% loss.

Of course I would need this coverage not directly over the coastline but land inwards. Is there any way to calculate this in ranges of percentage for the whole area I am studying?

The only way I see is to create grid of point data manually (using the same steps as here Calculating relative share of area covered by water/percentage of polygon area that overlaps in ArcGIS for Desktop?) and then making an interpolation of the data.


texts: Is it possible to display labels in background files (.imp) on Trimble GPS units?

texts:
I have been trying to figure out if this is possible for a while now. I know if you load a file in as a data file, you can click on a feature and the label will pop up. However, you cannot view a data file if you are collecting data and creating a new data file.
Has anyone had any success making labels show or know that it's not possible?


texts: How to Save as in 1 GeoTIFF file the hillshade raster and colored raster layer in QGIS 3.4.4?

texts:
I have a raster GeoTIFF of western Mediterranean sea and i want to apply the hillshade model with colors as the tutorial explains (Create great looking hillshaded maps in QGIS). After i finished my tests with hillshade and colors i want to save as 1 raster file both of layers ?

Can you help me ?


texts: looking for open source Java SE map viewer library that work with osm files or osm pre-downloaded tiles

texts:
I am new here and I hope I can find answers to my question here! I am developing a mapping platform for desktop, web, mobile (Android). I finally finished embedding mapping into my android software using osmdroid library. Now I want to start Java desktop part and I want to know is there any open source library like osmdroid or mapsforge for Java SE?

Thanks in advance.


texts: Layer Stacking landsat 5 surface reflectance high level data

texts:
I downloaded the surface reflectance landsat 5 scenes but when I layer stacked them my image was not how I imagined. I manually used atcor on normal scenes and it looks nothing like the scenes I layer stacked. Is there a certain order the tiffs have to be stacked for my image to come out correct?

I recently deleted the images and am going to try to use atcor myself. Do I need to download a calibration file for Landsat 5TM or can I use the default tm_essen.cal?


texts: I have a raster layer in ArcGIS with specific extent. How can I be able to find the angle of this raster frame? Because I create fishnet but it is like a rectangular area and does not match to direction of each pixel of this raster layer. So, I need correct angle to rotate it.


texts: How to determine if an image file is georeferenced

texts:
I have a task to scan a folder for georeferenced images. There might be a lot of images, some quite large, and some not georeferenced. The spatial information can also be either embedded or in a world file.

Programmatically (C#/WPF/ESRI Runtime) how can I tell if "C:\someFolder\file.x" is georeferenced.

Thanks


texts:
I have a very large raster stack of 6253 raster layers and I want to reclassify it.
For the reclass, I have tried the following but with no success:

my_rules <- c(-0.500, -0.251, 1,  -0.251, -0.101, 2,  -0.101,0.99,3, 0.99,0.269,4, 0.269,0.439,5, 0.439,0.659,6, 0.659,1.300,7, 1.300,2,NA)
reclmat <- matrix(my_rules, ncol=3, byrow=TRUE)
my_stack<- reclassify(my_stack, reclmat,right=FALSE)

Error in .local(x, rcl, ...) : all(rcl[, 2] >= rcl[, 1]) is not TRUE

Is this the wrong way to approach this?


texts:
I'm trying to implement a couple of Python actions so that you can start editing an object and save the changes from the id form.

I have had no problems adding the start edit action. However, saving the changes does not save the modifications that are made to the form.
from qgis.PyQt import QtWidgets
from qgis.gui import QgsAttributeForm

project = QgsProject.instance()
vl = QgsProject.instance().mapLayer('[% @layer_id %]')

vl.commitChanges()

QtWidgets.QMessageBox.information(None,u'Edición de registro [% $id %]'
,u'Registro guardado')

I have looked for the option to retrieve the form values but I have not succeeded.
Some help?


texts: How do I import OS raster tiles into QGIS at the correct scale?

texts:
I have some tiles of OS raster mapping and their associated TFW georeferencing files. The GR files are in the same folder as the raster data and have the same names, ie TQ21.tif/TQ21.tfw.

In QGIS 1.8, I add a tile as a raster layer and choose the OSGB 1936/British National Grid CRS.  When I load the next tile over, it loads in the correct place. Fine so far.

However, the problem is that the scale is completely out - to view a tile I have the onscreen scale indicator at around 1:365000000! If I go to Raster > Miscellaneous > Information there is nothing in the Raster Info box, but the fact that the tiles load in their correct place suggests QGIS is reading the GR files.

Am I missing something?


texts: What strategies, criteria, or rules to use for selecting coordinate systems?

texts:
What strategies, criteria, or rules do you use for selecting coordinate systems for

(a) storing,
(b) analyzing, and
(c) displaying GIS data?

(I humbly offer my reply to a related question about watershed analysis as an example of the considerations involved in (b).)

What are the pitfalls to watch out for?

Links to Web sites you find particularly helpful in this regard are welcome.


texts: Changing from a vector to a clear line (shoreline mapping)

texts:
I made a classification of land and water on an island. Especially on inner bays it becomes really confusing as the water is shallow. Is there a way to trace  line that shows the most likely shoreline (possibly with some kind of threshold)? To clarify I have attached one of the images.
I am using QGIS.


texts: How to retrieve the year from a date format?

texts:
I'm attempting to retrieve the year from the following date format: 1/1/1988
I've tried the following but it does not work str(1/1/1988)[4:] and I've also tried str(1//1//1988)[4:]


texts: Updating linestring geometry by replacing its points

texts: I have a linestring on my postgis. and I want to update a point in this linestring. I mean I want to update the geometry of the point in the line string by an other points


texts: How to export projected layers with relative cartesian (x,y) coordinates in ArcGIS

texts:
I've got several GIS datasets for a national park in South Africa (elevation map, water points, elephant density and such). I have access to the full ArcGIS Desktop Advanced suite, but am no GIS expert at all.

I want to use this data as a source of information inside of a multi-agent simulation. Therefore I need the coordinates of the resulting map in ArcMap to be transformed to relative x,y coordinates with a 0,0 origin somewhere instead of degrees and arc seconds (because programming and visualizing with those coordinate sytems is not really intuitive to me and the tools I use).

Actually I am ok with dealing with different file formats like ShapeFiles and .asc raster data in my simulation (using SharpMap), but I need the coordinates of these files to be relative, cartesian coordinates.

Can somebody tell me how to do that? I already got lost in the ArcGIS documentation.


texts: Arcpy: what tools are most suited for this job?

texts:
I'll try to be as much detailed as possible but please let me know if anything is missing:

I'm trying to get the number of lines in layerA that intersect with each line in LayerB and stores that number of intersections in a field in LayerB for each line, or row in tabular format. Both layers are huge and complicated loops may be problematic, but share your opinion please

the following is my logic:

-For each row in LayerB's Attribute table :
->Select(currentRow)->makefeature(currentRow)->selectByLocation(all lines that intersect with currentRow)->getCount(selected)->CalculateField(getCount)

will that theoretically work, please help me as I'm new to both arcGIS and python.


texts: Creating point feature class from multiple points using list comprehension and ArcPy?

texts:
import arcpy
sr = arcpy.SpatialReference(4326)

# Create polygon 2 x 2 with LL corner at 0,0...
coords = [[0, 0], [0, 2], [2, 2], [2, 0], [0, 0]]
polygon = arcpy.Polygon(arcpy.Array([arcpy.Point(*a) for a in coords]), sr)
arcpy.CopyFeatures_management(polygon, r"C:\Temp\polygon.shp")

# Create point inside the polygon...
point = arcpy.PointGeometry(arcpy.Point(1, 1), sr)
arcpy.CopyFeatures_management(point, r"C:\Temp\point.shp")

I can quickly create a polygon object using a nested list comprehension statement within the arcpy.Array function and use it as a parameter in the CopyFeatures geoprocessing tool (as shown above).

I can quickly create a feature class with a single point object also shown above, again using the CopyFeatures geoprocessing tool.

Is it possible to create a "point" feature class from the "coords" list variable using a similar nested list comprehension approach?

For example, I tried with the following but obviously something is not right here.

point = arcpy.PointGeometry(arcpy.Array([arcpy.Point(*a) for a in coords[:-1]]), sr)

Just trying to make a quick one-liner without resorting to the arcpy.da.InsertCursor method if possible.


texts: ArcGIS Geometric Network splits at junctions after editing route?

texts:
When I edit a Vertex of route (1) in my geometric Network, the route (2) which first was one route splits up in two routes at the junction. The m-Values from the Initial route are from 0.000 - 255.345. After the split the values are 0.000 - 145.673 and 145.673 to 255.345.
Anyone knows why the route gets split up in a geometric Network? I thought the Network should ensure connectivity?
My Network is build correctly.


texts: Converting from local coordinates to coordinates GPS?

texts:
i am working with business object tool traying to make a kml file to be imported on google maps. I dont know much about geographic systems, but i begin to understand some basic things.

I have  a database which give me a local coordinates i know that because the gps has an appfile which has a grid configured and the database says me that units are meters.

The app file shows this information:

Site Point Pairs:
        Number Of Points:   None

Ellipsoid:
    Name:   (WGS 84)
    Semi Major Axis:    6,378,137.000
    Semi Minor Axis:    6,356,752.314
    Flattening: 1/298.2572215381
    Eccentricity:
        First:  0.0818191911
        Second: 0.0820944382

Datum Transformation:
    Name:   (WGS 84)
    Method: Molodensky
    Parameters:
        Direction:  WGS84 to Local
        Translation X:  0.000
        Translation Y:  0.000
        Translation Z:  0.000

Geoid Model:
    Name:   Subgrid of MEXICO97 (Mexico)
    Method: Grid
    Geoid Grid File Name:

Coordinate System:
    Name:   12 North
    Group:  UTM

Projection:
    Type:   Mercator Transversa
    Parameters:
        Latitud central:    0°00'00.00000"North
        Longitud central:   111°00'00.00000"West
        Falso norte:    0.000
        Falso este: 500,000.000
        Factor de escala:   0.9996
    AZI Parameters:
        Azimuth Direction:  North Azimuth
        Positive Coordinate Direction:  NorthEast
Embedded Data:  True
Embedded Geoid Grid Extents
Minimum Latitude:  30°53'60.0000"North
Maximum Latitude:  31°18'0.0000"North
Minimum Longitude: 247°03'60.0000"East
Maximum Longitude: 247°33'60.0000"East

CSIB Version: 10.50
CSIB Format: Big Endian

the database shows me this kind of values

vehicle     x          y          z

VE-28      328279     3457947    -48
VE-28      328295     3459609    192
VE-28      328288     3459617    193
VE-28      328304     3459622    192
VE-28      328220     3457362    -112

I would like to convert this unit into gps information through some formula.
or get the gps to give the gps coordinate to the database instead of the local coordinate.


texts: How do I reassign a point's XY value in ArcMap 9.3?

texts:
I have a rectangle polygon I drew on a map and then added points to it's corner and midpoint. I drew the shape (representing a forest sample plot) where I thought it was best located based on 70 year old hand drawn maps. However, I now know the exact GPS coordinates (UTM) of the plot and would like to move the shape and points on my map to the correct GPS points.

Is there a way to select these features and change their XY coordinates to match the UTM values I now have?


texts: How to share Leaflet web map with local data

texts: I create a web map using Leaflet. I access my map through localhost. However now I want to share with other persons. But I have local data. If I open the HTML it only presents the base maps. How can I share my local data?


texts: Moving rasters from ArcGIS 10.0 to Mapinfo Pro 15.0?

texts:
I have several rasters in ArcGIS in GRID format. I need to move them into Mapinfo environment and I'm a bit confused about the most optimal format. Overall, is it recommended to export ArcGIS GRID files into some other format when moving to Mapinfo? And what might the best one be? Rasters will be used just as a background, no further analysis needed.

For example when opening these rasters in .TIF format, all I get is a black image following image extension.


texts: How to get OSM data for different areas at once?

texts:
I am trying to get some OSM data (schools, shops, public transport stops etc.) for several cities. I have a shp-boundaries for cities.

I use QuickOSM for QGIS or overpass-turbo on a Web. I have a query for data but I can't figure out how to automatize extraction for different bounding boxes at once. For now I select bounding box every time manually :(


texts:
I have a full planet imported in PostgreSQL with imposm3.

What is the best way to generate vector tiles (MVT in .mbtiles) from PBF/PostgreSQL with maximum hardware utilization (distributed generation to multiple CPUs/cores)?


texts: Attempting to map GeoJson file producing issues

texts:
I've been trying to get this to work for some time now. The end goal is to produce a map I can show using d3.js, however d3 is butchering the map and producing an undesirable image. The code is as follows:

var width  = 1000;
var height = 800;

var vis = d3.select("body").append("svg")
  .attr("width", width).attr("height", height)

d3.json("ca-ns-all.geo.json", function(json) {
// create a first guess for the projection
  var center = d3.geo.centroid(json)
  var scale  = 150;
  var offset = [width/2, height/2];
  var projection = d3.geo.mercator().scale(scale).center(center)
      .translate(offset);

  // create the path
  var path = d3.geo.path().projection(projection);

  // using the path determine the bounds of the current map and use
  // these to determine better values for the scale and translation
  var bounds  = path.bounds(json);
  var hscale  = scale*width  / (bounds[1][0] - bounds[0][0]);
  var vscale  = scale*height / (bounds[1][1] - bounds[0][1]);
  var scale   = (hscale < vscale) ? hscale : vscale;
  var offset  = [width - (bounds[0][0] + bounds[1][0])/2,
                    height - (bounds[0][1] + bounds[1][1])/2];

  // new projection
  projection = d3.geo.mercator().center(center)
    .scale(scale).translate(offset);
  path = path.projection(projection);

  // add a rectangle to see the bound of the svg
  vis.append("rect").attr('width', width).attr('height', height)
    .style('stroke', 'black').style('fill', 'none');

  vis.selectAll("path").data(json.features).enter().append("path")
    .attr("d", path)
    .style("fill", "red")
    .style("stroke-width", "1")
    .style("stroke", "black")
});

And the map is defined from the data below

http://code.highcharts.com/mapdata/countries/ca/ca-ns-all.geo.json