texts:
I have a big round polygon feature (5km diameter) and a small round polygon feature (3km diameter) that I would like to clip with Python GDAL/OGR to obtain the "donut". Once this is done, I will apply the same process for 25.000 features, that is why I need to do this with GDAL/OGR and not editing manually in QGIS/ArcGIS. It seems that the Clip operation works at the level of layers. Therefore, what I did was creating a layer for my 5km circle, a layer for my 3km circle and an output layer required by the prototype of the function.

Apparently the operation is finishing correctly, because I get a .shp file with 6KB inside, but when I open this shapefile either with QGIS or ArcGIS it is completely blank. I thought it was the projection, but the .prj file says it is set to EPSG:28992 (Amersfoort_RD_New), which is what I need.

The code I wrote is the following:

outDriver = ogr.GetDriverByName("ESRI Shapefile")
outDataSource = outDriver.CreateDataSource(pathout)
outLayer = outDataSource.CreateLayer("Big", srs = outSpatialRef, geom_type=ogr.wkbPolygon)

outDriver2 = ogr.GetDriverByName("ESRI Shapefile")
outDataSource2 = outDriver2.CreateDataSource(pathout)
outLayer2 = outDataSource2.CreateLayer("Small", srs = outSpatialRef, geom_type=ogr.wkbPolygon)

outDriver3 = ogr.GetDriverByName("ESRI Shapefile")
outDataSource3 = outDriver3.CreateDataSource(pathout)
outLayer3 = outDataSource3.CreateLayer("Clip", srs = outSpatialRef, geom_type=ogr.wkbPolygon)

outLayer.CreateFeature(mybigfeature)
outLayer2.CreateFeature(mysmallfeature)
something = outLayer.Clip(outLayer2, outLayer3)
print something # returns zero

So, I do not understand how the operation seems to finish correctly, writes 6KB in my disk and then I see nothing on QGIS. I have attached an image, showing the features, but not the clipped polygon.

I guess there is something missing: memory flush?

insert a geometry for the output (how if there is no return of Clip)?

close the files?

Hope someone can spot the problem!

Following the advice of Benjamin, I installed Shapely and Fiona libraries to solve this problem. Shapely is used "for manipulation and analysis of planar geometric objects" and Fiona does the reading and writing of shapefiles in a Pythonic and very neat way.

So after playing around with these libraries, I came up with the following code:

from shapely.geometry import shape, mapping
import fiona

pathb = r'PATH_TO_BIG_ROUND_FEATURES\Buffer_5km.shp'
paths = r'PATH_TO_SMALL_ROUND_FEATURES\Buffer_2km.shp'
pathout = r'PATH_TO_OUTPUT_FOLDER\clip.shp'

big = fiona.open(pathb)
small = fiona.open(paths)

polbig = big.next()
polsmall = small.next()

a = shape(polbig['geometry'])
b = shape(polsmall['geometry'])
c = a.difference(b)

source_schema = { 'geometry': 'Polygon', 'properties': {'id': 'int'}}
source_crs = {'no_defs': True, 'ellps': 'WGS84', 'datum': 'WGS84', 'proj': 'longlat'}
source_driver = "ESRI Shapefile"

with fiona.open(pathout, 'w', 'ESRI Shapefile', schema=source_schema, crs=source_crs) as out:
    out.write({
        'geometry': mapping(c),
        'properties': {'id': 123},
    })

And this produces the following output:

Which is exactly what I need!


texts: How to Blink a line on a map similarly like a point blinking on map through lat lon points

texts:
how to blink a line on a map which is formed by lat lon points through linestring from database the points(lat,lon) comes which is stored in form of MultiLineString.
similarly like how a point blinking on a map same way the line should blink on map.
(or)

I need to show on a map line blinking below code is displaying a line on the map but it should be made line blinking not all but particular some lines.

code :-

var mainfed = new OpenLayers.Layer.Vector("Mainfed");

$.get('http://'localhost:90'/them/date-feeder.php',function(data) { //getting the data from DB through php

   var records = data.split(',<br>'); //recognize the lat,lon through splitting the  records

 for(var idx=0;idx<records.length;idx++)
  {
    var myarr = records[idx].split(",");
    var points = new Array()
    for(i = 0; i < myarr.length; i++)
    {
        var lon = parseFloat(myarr[i++]);
        var lat = parseFloat(myarr[i]);

                      points.push(new OpenLayers.Geometry.Point(lon,lat).transform(new OpenLayers.Projection("EPSG:4326"), map.getProjectionObject()) );
    }
         var line = new OpenLayers.Geometry.LineString(points);
         var style = {
         'class': 'lineblink',
          strokeColor: 'red',
          strokeOpacity: 5,
          strokeWidth: 2
          };
 var lineFeature = new OpenLayers.Feature.Vector(line, null, style);
    mainfed.addFeatures([lineFeature]); // forming a line through points
         }
});
    map.addLayer(mainfed);
map.addControl(new OpenLayers.Control.DrawFeature(mainfed, OpenLayers.Handler.Path));


texts: Merging seemingly contiguous polygons leaves artifacts

texts:
I am working with a NaturalEarth country vector dataset (download it here) and I am trying to merge country polygons according to custom regions (based on the UN Geoscheme), in order to only have an outer border rather than every single country's border.

My issue is that, seemingly contiguous countries leave artefacts after merging, be it by hand ("edit > merge selected features") or by using the "dissolve" tool ("vector > geoprocessing tools > dissolve").

How come there are stray lines inside the regions?
Can I somehow fix that, or should I find a different dataset?

See the following examples in Africa:

Merging by hand

Dissolving according to data


texts: WKT for EPSG:4326 with LON 0 to +360 instead of -180 to +180

texts:
I have GRIB2 files that seem to be in a projection similar to EPSG: 4326, except that instead of longitude between -180 and +180 it is between 0 and 360.

My intention is to get these files into geoserver and allow the data to be queried and reprojected at will through the WMS. I can get these files into geoserver with the GRIB extension, but I get a null pointer error when I attempt to reproject to mercator. I assume that this error occurs because data lies outside the -180 to +180 longitude bound.

I had attempted to just reproject the GRIB files with gdal_translate, using the PROJ4 +lon_wrap=180option, but that operation takes forever and produces a far larger output file for some reason.

From what I understand, I can fix this situation by creating a custom projection in user_projections/epsg.propertiesthat describes the CRS of these GRIB files. A custom projection for geoserver must be in the WKT format.

TL;DR: How do I write a WKT that describes a lat/lon coordinate system similar to EPSG: 4326 but with longitude ranging from 0 to +360 instead of -180 to +180?

Here is the gdalinfo for one of these GRIB files. The file depicts the continental United States:

Size is 2847, 1517
Coordinate System is:
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (229.982699929725896,55.001944129287601)
Pixel Size = (0.024600140548138,-0.023088258575198)
Corner Coordinates:
Upper Left  (     229.983,      55.002) (229d58'57.72"E, 55d 0' 7.00"N)
Lower Left  (     229.983,      19.977) (229d58'57.72"E, 19d58'37.40"N)
Upper Right (     300.019,      55.002) (300d 1' 9.48"E, 55d 0' 7.00"N)
Lower Right (     300.019,      19.977) (300d 1' 9.48"E, 19d58'37.40"N)
Center      (     265.001,      37.490) (265d 0' 3.60"E, 37d29'22.20"N)

gdalsrsinfo:

PROJ.4 : '+proj=longlat +a=6371229 +b=6371229 +no_defs '

OGC WKT :
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]

Screenshot from QGIS:

Edit: This is not a duplicate. I have seen the answer from the referenced question, and it is not satisfactory. The operation takes too long to complete and produces a file that is entirely too large.

Edit #2: Added a link to the GRIB file:
 http://furlender.com/ECMWF_2017052412.grib2


texts: What is fix geometry actually doing in QGIS?

texts:
Sometimes geometry's have to be fixed for further processing. So you can use fix geometry and it will do the job for you without prompting. But I'm wondering, what the tool fix geometry is really doing.
What kind of geometry errors is it fixing?


texts: Is the layer value of MODIS Day_view_time and Night_view_time in Local Time or in Local Solar Time?

texts:
On the NASA homepage in the layer description you can find following line for the Day_view_time layer:

Does this mean that the layer contains the Local Time (LT) in the corresponding local time zone or in Local Solar Time (LST)?
And if this is the LST, is the difference significant if I'm not interested in the exact minute of the observation time?


texts: Adding spatial reference to many shapefiles in loop using ArcPy?

texts:
AST folder has 5 shapefiles.
This code tries to pass the 2100 coordinate system to these shapefiles. But it gives error 999999:

arcpy.env.workspace = r"C:\Users\user\Desktop\roads\AST"
for i in arcpy.ListFeatureClasses():
     sr = arcpy.SpatialReference("EPSG:2100")
     arcpy.DefineProjection_management(i, sr)

How to make it work?


texts: Is there a way to retrieve all cities in the world from OpenStreetMap?

texts:
I have tried to retrieve them from the api with the following:

wget -O allcities.xml "http://www.overpass-api.de/api/xapi_meta?node[admin_level=8]"

however it does not seem to return all of them. From the doc Admin Level 8 should be cities, towns and villages but it seems to only contain major cities. I get a list of about 12000 nodes.

I am actually looking into getting all inhabited places on earth with at least 1000 population.

I am also downloading the whole planet in PBF format, in case using the API is not possible. How would I be able to achieve the same using such file?


texts: R Raster layer of pixels overlapped by shapefile of points

texts:
Given a raster layer and a shapefile of points.
How can I get, in R, a new raster layer of pixel values only where there are the point, while all other raster pixel values are kept as NaN?


texts: WCS : InvalidParameterValue and OperationNotSupported : Differences

texts:
In WCS 2.0 & 1.1, I'm looking for generate exceptionCode "OperationNotSupported" which is defines in OGC Documentations (WCS 1.1 / WCS 2.0 / OWS 1.1 / OWS 2.0). With MapServer 7, I've tested: http://hostname:port/path?service=WCS&request=GetZZZ.

XML response gives en error with exccode "InvalidParameterValue"
I understand why there is this exceptionCode because it's an invalid parameter but when I can do the output with exceptionCode "OperationNotSupported"?

What do I do for that?

What is the difference between "OperationNotSupported" and "InvalidParameterValue"?


texts: Open dataset(s) with centroids or other geometry for WOEIDs?

texts:
The Yahoo! GeoPlanet data is available for download but this data only includes the object hierarchies without any associated geometry.

Do anyone know if there is an open source for geometry (centroids or otherwise) for associated WOEIDs? I wonder if there's a lesser-known link somewhere on the GeoPlanet site, or perhaps another party has done some work to correlate open data sets (census, geonames, etc) with WOEIDs.

A WOEID (Where on Earth IDentifier) is a unique reference identifier assigned by Yahoo! to identify any feature on Earth.


texts: Installing xlutils alongside xlrd and xlwt in ArcPy install folder to facilitate copying Excel worksheets?

texts:
Whenever possible I avoid installing Python libraries other than those installed with ArcPy by the ArcGIS Desktop 10.x and ArcGIS Pro installs but occasionally I have for ReportLab, gdal, pyodbc and I usually find several ways are suggested and I only get them working after some experimentation.

I would like to be able to copy a worksheet from one Excel (*.xls) spreadsheet, created using TableToExcel in ArcPy, to another Excel (*.xls) spreadsheet (also created by ArcPy), and the only technique that I found that looks promising for doing this is in this answer to Writing to existing workbook using xlwt

That technique commences by importing xlrd, xlwt and xlutils using the re-ordered lines below:

from xlrd import open_workbook
from xlwt import easyxf
from xlutils.copy import copy

When I run them using IDLE from Python 2.7.10 that was installed with ArcGIS 10.4.1 for Desktop I get this error which indicates that xlrd and xlwt are available but xutils is not.

  Traceback (most recent call last):   File "C:\Fyfe\xutilsTest.py",
  line 3, in
      from xlutils.copy import copy ImportError: No module named xlutils.copy

This matches what I see when I look in C:\Python27\ArcGIS10.4\Lib\site-packages where xlrd and xlwt folders can be seen but no xlutils.

Is there a simple process that will let me download a folder of Python code for xlutils and copy it into the same folder to make import xlutils just work?

Installing third party python libraries in ArcGIS location gives me some hints about what may need to be done but that Q&A seems to be more understandable by someone who develops using a higher Python:ArcPy ratio than I do.


texts: Using letter and word spacing unit in QGIS labeling?

texts:
I just want to know if letter and word spacing, in QGIS labeling, are expressed in points, and if the value its absolute or relative to the font.

I have searched a lot, but I can't find any clarification about it.


texts: Why does the Service information disappear on ArcGIS Server Manager?

texts:
I have 41 services on our ArcGIS Server. The Server Manager holds roughly 10/page. When I click on any page other than page 1, the service information disappears (status, number of instances, etc.)

When I refresh, it always takes me back to page 1. Any ideas on how to fix this issue?

Running ArcGIS Server 10.1

The services disappear in IE 9 and Firefox. I can't login via Chrome so no way of telling if it would work or not.


texts: How to make a Python script that does a task in QGIS?

texts:
When I run this code in QGIS Python console it works fine, but I want want to run it from an external Python script.

What I want the script to do is:
1- Open QGIS
2- Run this code
3- Close QGIS

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtXml import *

mapinstance = QgsMapLayerRegistry.instance()

uri =  "file:///C:/Users/Home/Desktop//sample.csv?delimiter=%s&xField=%s&yField=%s" % (",", "Long", "Lat")

vlayer=QgsVectorLayer(uri,"sample","delimitedtext")
mapinstance.addMapLayer(vlayer)

uri="C:/Users/Home/Desktop/Style.qml"
vlayer.loadNamedStyle(uri)

vlayer.selectAll()
qgis.utils.iface.actionZoomToSelected().trigger()
vlayer.invertSelection ()

qgis.utils.iface.mapCanvas().saveAsImage('/Users/HOME/Desktop/capture.png')


texts: Where does ArcGIS Desktop get image data from?

texts:
How does ArcGIS Desktop (and probably most other WMS clients) request data from supported servers? If we take their own sample server as an example, I can understand how Desktop connects to the server, knows which layers are on it etc, but when I add the service to Desktop and it starts showing map images, how does it get those images? I can't find anything in the XML file that points to any kind of source for these images, although it has SRS, boundingbox and all the metadata about the maps.

What am I missing here?

EDIT: To clarify, I have a (homemade) server application, which creates and sends the following XML file upon receiving a GetCapabilities request:

<?xml version="1.0" encoding="UTF-8"?>
<WMS_Capabilities xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
    <Service>
        <Name>Name</Name>
        <Title>EivaTMS</Title>
        <Abstract>Something clever about this WMS</Abstract>
        <Capability>
            <Layer>
                <Title>kraken.db</Title>
                <Name>Name</Name>
                <Abstract/>
                <SRS>OSGEO:41001</SRS>
                <Profile>NaviPacLocal</Profile>
                <Format>image/png</Format>
                <BoundingBox>
                    <xmin>140</xmin>
                    <xmax>140</xmax>
                    <ymin>80</ymin>
                    <ymax>80</ymax>
                    <zmin>0</zmin>
                    <zmax>0</zmax>
                </BoundingBox>
                <MapSize>256</MapSize>
                <CellSize>3</CellSize>
                <MaxLevel>32</MaxLevel>
            </Layer>
        </Capability>
    </Service>
</WMS_Capabilities>

I'm well away that this XML doesn't stand a chance when comparing it to the WMS specification but this is due to it being used for my own client at some point, but now I'm in need of being able to use it with ArcGIS. Now, ArcGIS Explorer is obviously able to send a GetCapabilities request (some_url?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities) and it's able to see 'kraken.db' which I assume is because of the 'Layer' tag it's surrounded by. I can see my server receives the request and that it sends the XML back to ArcGIS, which gives me this window where I have to choose a service:

But when I try to add it, I get a "Unable to add service(s)" error. This isn't strange since my XML lacks a lot of information, but I have difficulties finding out exactly what kind of information ArcGIS Explorer needs in this situation.


texts:
I am running Ubuntu 18.04 with the XFCE interface and QGIS v3.6. When I start QGIS the menu bar at the top of the window is missing, and I can't find any way to bring it back. I found some solutions online that suggest modifying the QGIS2.conf file, but I don't have one of those.

Do any of you have ideas how to solve this? The main reason I use QGIS is to save quick previews of images, which is no longer possible without the menu bar!

(edit) Rostranimin suggests launching from the command line with the no configuration option, I had actually already tried this (but neglected to mention it). Doing so still results in no menu bar.
(edit2) To give more background info:
I was previously using Ubuntu 14.04 on this machine and whatever the default version of QGIS was for that distro (2.4?), then my company upgraded me to 18.04 with QGIS 3.6.
In 2.4 the menubar was visible, it hasn't been visible since the update. Unfortunately as this is a company-managed machine I don't have sudo access, so my options to change things are limited.

I have another machine (self-managed) that's running Mint 19 with QGIS 3.6, the menu bar appears fine there.


texts: Summing up the area on a polygon by a certain attribute result

texts:
I wish to know if there's a tool on QGIS equivalent to the Excel's "SUMIF".

I have a shapefile of a land use and cover that's classified by the attribute with the name of the classes, for example: "vegetation", "urban centers" and "pastures". My goal is to calculate the area of each class, but since it's composed by various polygons, I need to add up all the areas from all polygons that were put in the same class.

Like I mentioned in the beginning of the question, I know that Excel has a function that can do such thing for me, but it turns out that the sum made by excel doesn't match the area's total given by QGIS.

English is not my first language, so if the explanation isn't clear enough, I'll be glad to rephrase it.


texts: Combining multiple .tiff files obtained after download from Google Earth Engine?

texts:
I downloaded some of my classified images from GEE using the following script-

Export.image.toDrive({
  image: input.clip(district),
  description: 'Test',
  maxPixels: 499295920080,
  scale: 30
});

But since my image was large, it was broken down into 16 .tiff files during download. Now I am required to do some post processing on my complete classified image as a whole, but I can't figure out how to combine these images back into a single jpeg/.tiff image.

Can someone guide me towards a viable solution?


texts: How do I know when to use a map server or a tile server?

texts:
I'm new in GIS tech and I still have some questions about GIS.

I wanted to know why to choose a map server (MapServer, geomap) instead of a tile server (tilestream, tilecache or tilestache)  + tile renderer (Mapnik)?

I know that a server map use spatial data from a DB, it generate image that can be read by client map library (OpenLayers...).
For example if I have a shapefiles and I don't need dynamic data, what the point to use a server map?

There is no need if I have a static map to have a server map, is that true?

When does a map server regenerate map, does it have a cache?

What is the format of the generated tile ?

If I choose a map server, how can I style my map? does it have a soft like in Tilemill or do I have to do it directly in the code?

Does anyone know software which does coloring like tilemill? or do I have to use cartoCSS or cascadenik to do it?

I a bit confused about that, so if anyone got the answer or think I misunderstood something, it would be great to share it!


texts: How to create "roof line" for a general polygon area?

texts:
I have a polygon with an arbitrary shape ( assume that it is well-formed with no self intersection, no duplicated points etc), and I need to subdivide the polygon into sub-polygons like how one would define roof on top of a building.
How can I do that in any GIS software, or through any open source libraries, or through any algorithm?

The reason why I'm asking this is because for my flat terrain ( with 2D shape unknown obviously), I have to manually enforce artificial gradients at appropriate places, so that the water can flow to the neighboring cells correctly when I am applying D8 implementation to calculate the flow accumulation. So the most intuitive places for sub-division would be the places where roof lines are taking along.

I would be interested to know is there any algorithm that one can use for this, and also is there any C++ open source libraries implementation for this ( I don't want to reinvent the wheel); is GDAL or CGAL or any geometrical algorithm already has this implementation?


texts: What's the best way to get the OBJECTID name?

texts:
What's the best way to get the OBJECTID name of a feature class?

I know you can use "OID@" in a cursor to get the field, but can you get the name from that as well?

I'm writing a script, and need to take into account varying OBJECTID names, like OBJECTID_2, OBJECTID_12, FID, etc. that will be used in a query.

Is using Describe the only way to get the name?

For example:

oid_field = [field.name for field in arcpy.ListFields("FC") if field.type == "OID"]

I'm just wondering if there's another way.

Thanks.


texts: ASCII to RASTER failing with comma syntax error?

texts:
ASCII to RASTER fails with comma syntax error, but data appears correctly formatted with no commas

I'm trying to use the ASCII to RASTER tool (ArcMap V10.4.1) on PRISM rainfall data.  The file appears to have the correct format for this use:

ncols 1405

nrows 621

xllcorner -125.020833333333

yllcorner 24.062500000000

cellsize .04166666666667

NODATA_value -9999

-9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00

... and so on...

However, I get the following error messages:

ERROR 010328: Syntax error at or near symbol ,.

ERROR 010267: Syntax error in parsing grid expression.

Failed to execute (ASCIIToRASTER).

There are no commas in this entire file.  So I'm at a loss for what this error refers to.

Full PRISM asc files can be downloaded here: http://www.prism.oregonstate.edu/recent/


texts: proj4text column blank for SRID 97688 - causing errors

texts:
Apparently a blank proj4text field causes errors when adding a record. When you enter something close to correct into proj4text, the records get inserted without error.

What is the correct value to use, since the one I cobbled together gave bad results?

This is the Insert record from SpatialReference.org for that SRID:

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 97688, 'sr-org', 7688, '',
'PROJCS["WGS_1984_Web_Mercator",GEOGCS["GCS_WGS_1984_Major_Auxiliary_Sphere",DATUM["D_WGS_1984_Major_Auxiliary_Sphere",SPHEROID["WGS_1984_Major_Auxiliary_Sphere",6378137.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1.0]]');

What is the correct value for proj4text given that srtext?


texts: How to determine if all points on a map are identical in arcpy?

texts:
My eventual goal is to automate kernel density calculations for a large corpus of point data; however, if all points are identical (for example, if there are 7 points all stacked on top of one another because they have the same long/lat values), the kernel density function in arcpy will not work. It will result in the following error:

"Error 010246: All input points are identical A valid window must be specified"

Since this will be an automated process, I need a way to check the points in the current input to be sure they aren't all the same using python/arcpy. Is there a way to do that?

EDIT: I specified a grid size and an extent, so a lot of the errors are fixed; however, identical points still give me the "valid window must be specified." I have tried to use the FindIdentical function to look for identical points, but this returns a DBF file, that I cannot read in the script. Is there some way to read a DBF easily in python, or can I find if all points are identical in another way?


texts: How to remove, using script, non present entry inside legend symbology?

texts: Using Qgis after loading shp file i loaded a general qml files that are coded for a specific field for my shp.  Inside qml are coded a larger number of symbology of the field that not all are present inside my shp. I will from pyhton console delete all entry in legend (loaded from qml) that have zero count geometry inside the shp, suppose Shp to purge is the current.


texts: Geowebcache alias for custom projections

texts:
Is there a way to configure a projection alias in geowebcache?

E.g. if I would like to advertise projection EPSG:25833 and EPSG:32633 to the client, while requesting tiles from the backend using EPSG:25833. That is for (mostly) identical projections where no reprojection is needed.

The idea being of course that identical tiles are only cached (and requested) once, while both projections are advertised.


texts: Why is GPS almost unavailable without data connection?

texts:
When travelling to foreign countries, ie. when a data connection is unavailable because I don't want to spend big bucks using a roaming data connection, I noticed that GPS is virtually never available.

It can take hours for GPS to show my location, and even then, only when I'm located in open areas for a long time, with no buildings or trees. The only time GPS is available is when the phone is connected to the Net through wifi, which is pretty much useless.

I am using a Galaxy Nexus running CyanogenMod, and the phone has the FasterGPS utility installed.

Can GPS be available in cities without a data connection?


texts: My address is not shown in the maps how to update the maps with my new address

texts: I have moved to my new house and not able to find it on Google maps, how to add my address to the maps.


texts:
I did not realize that I should have been deleting cursors after accessing attribute tables until today. Up until now I never did this, but I've noticed that there are lots of file locks after my script runs (I presume it was a bad practice on my part to begin with). ArcGIS 10 help tells me that I should use the following code:

# Delete cursor and row objects to remove locks on the data
#
del row
del rows

I'm doing that in the following function:

# calculate areas
def area_calc( shapefile ):
    "calculates geometric properties of each building"

    # loop through attribute table
    Rows = ARCPY.UpdateCursor( shapefile )

    for row in Rows:

        geom = row.Shape

        polyArea = geom.area
        row.setValue( perimName, 0 )  # set to zero
        polyPerim = geom.length

        row.setValue( areaName, polyArea )   # update area value
        row.setValue( perimName, polyPerim ) # update perimeter value

        Rows.updateRow(row)

    # Delete the row and cursor
    del row
    del Rows

    ARCPY.AddMessage( 'Area and Perimeter calculated for each polygon' )

However, when I do that I get the following error message:

Traceback (most recent call last):
  File "LoopModuleBuildingAnalysis.py", line 9, in <module>
    import ConstructionMaterialPhase1
  File "D:\Code\AnalysisPhase1.py", line 101
    del row
      ^
SyntaxError: invalid syntax
Press any key to continue . . .

Any suggestions what I should be doing?


texts: How to create customized side control bar in OpenLayers?

texts:
I could create a simple map using OpenLayers. But I want to have a customised side control bar.

The one I created:

The one I want to create:


texts:
I have tried performing the query shown below via ogr2ogr and the gpd.GeoDataFrame.from_postgis() method in geopandas.

from osgeo import ogr
import psycopg2, os

import geopandas as gpd
import pandas as pd
ogr.UseExceptions()

server = "xxxx.com"
name = "xxxx"
user = "xxxx"
pw = "xxxx"
port='5432'
lyr_name = 'ships'

conn = psycopg2.connect('connection details')
cursor = conn.cursor()
print('executing sql')
query = 'select * from ships a, ships b \
where a.dt_upload>\'2019-01-01\' and b.dt_upload>\'2019-01-01\' \
and a.det_id>b.det_id and ST_DWithin(a.ship_centroid, b.ship_centroid,100)'
#query executed below always fails
gdf = gpd.GeoDataFrame.from_postgis(query,conn, 'ship_center1', crs='EPSG:4326')

ships1=gdf.iloc[:,0:7]
ships2 = gdf.iloc[:,7:]
ships2['shipcenter2'] = ships2['shipcenter2'].apply(wkt.loads)
ships2.columns = ships1.columns
new_gdf = gpd.GeoDataFrame(pd.concat([ships1,ships2], ignore_index=True),crs = {'init': 'epsg:4326'} )
new_gdf.rename(columns={'shipcenter1':'geometry'}, inplace=True)
new_gdf.to_file('./ship_coverage.shp' )

Everytime I try the process gets killed. There are about a million rows in my table. How can I optimise this query to fully execute so I can either export the results as a shapefile or as a GeoDataFrame?


texts: ArcGIS Error: ERROR 000223: Cannot extract event table properties while trying to create an event layer using ArcPy

texts:
I have a .dbf file and want to change it to a .shp file with its LAT/LON.

ArcPy Code,

import arcpy

arcpy.env.overwriteOutput=True
dbffile = "GIS_counts.dbf"

spref =  arcpy.SpatialReference(4326) #gcs_wgs_1984
arcpy.MakeXYEventLayer_management (dbffile, "LONGITUD", "LATITUDE", "test.shp", spref)

Error:

ERROR 000223: Cannot extract event table properties
Failed to execute (MakeXYEventLayer).


texts: What is the difference between creating a plugin in QGIS 3.x and writing a python standalone script in QGIS 3.x?

texts: I am developing a plugin in QGIS where I have used the template "Tool button with dialog". I am using this tutorial: http://g-sherman.github.io/Qgis-Plugin-Builder/ and https://www.qgistutorials.com/en/docs/building_a_python_plugin.html . However, I came across this http://www.qgistutorials.com/fi_FI/docs/processing_python_scripts.html where a python standalone script in the processing toolbox will also easily create a similar user interface. Which one is more user friendly and what is the advantage?


texts: Still no fictional geographies in some GIS format?

texts:
I'm actually reviving Are any fictional worlds available in some GIS format? from 2011 as I'm hoping the Game of Thrones fever has boosted people's interest in drawing fictitious countries in some GIS format (e.g. .shp) from 2011 until now, so that others can improve or modify.

Alternatively, does anyone know of an easy way to transform a bitmap map such as the one below to a vector map?

Source: DeviantArt.


texts: Does shapefile or MapInfo TAB have advantage as source format for OGR conversion to GeoJSON?

texts: I'm relatively new to OGR and GIS in general. Most of the sources I look at give a choice of ESRI shapefile and MapInfo formats. I'm using ogr2ogr to convert these to GeoJSON for import into JSON-based databases. Can any experienced OGR users tell me if there is any advantage in the choice of a source format here? I'm more interested in preserving the integrity of the data than speed, etc.


texts: Sea level rise prediction by zip code in the U.S. coasts

texts: I'm relatively new to GIS data and am working on a project on climate change. I can find online predictions of sea level rise by zip code in the U.S. East, West and Gulf Coasts (for example:  http://sealevel.climatecentral.org/) but I would like to know if there is public data available so that I could generate similar predictions in python or a similar language.


texts: Measuring shortest distance along line feature between 2 polygons?

texts:
I am using ArcGIS and trying to construct the ModelBuilder to calculate the shortest distance along the line feature for a pair of polygon.

For example, as in the attached picture, I have polygon number 254 and 297 and the line feature (in red).

Is there any alternatives or commands in ArcGIS to help me measure the short distance on the red line traveling from polygon 254 to 297?

Edit:

The picture I send is just the simple problem to understand what I want. The real project is targeted on more than 700 polygon and full transport network of public transport.

So, the simple measure and editing manually is not possible.

to be a little clear what I have for the data is:
1. Public Transport Network
2. OD matrix table

what I need is identify daily trip passengers for throught trip (taking only one route of transit from O to D without transfering) for each transit route.

I did not consider to use network analysis because, first, I did not really familiar with it and skim through the network analysis tutorial, I couldn't find its ability to storing information of trip volume per route (or by ID).

what I have done so far is figure it out to using geo-processing script to automate the process of line selection/split/measure of the shortest length of the line between O and D.

what my model is now mostly look like the attached file and what I have a problem is now I am now able to select the route that traverse between to polygon (O and D) but I have to split and select the exact line segment for the shortest distance.

The model builder is  not allow to use select by attribute or location for the layer that not exist on the project windows.  So, the next process have to start all over again by adding the outputs of my current works to the next model builder program.

So, what I really want is just curious that do they have any simple tools in GIS to measure this simple problem in Model Builder???  distance on the line between two point or polygon boundary (****again not manually measure)


texts: Recalculate fields value using ArcPy and Calculate Field tool?

texts:
Im not sure how to pass through a list of fields within the calculate field tool.I have a list of fields and each of them need to be divided by 1000.
My code is as follows:

import arcpy
shp = r"m\team.shp"
fields = arcpy.ListFields(shp, "*qty")
for f in fields:
    print f.name
    arcpy.CalculateField_management(shp, f, "f/1000", "PYTHON")


texts: automatic zoom to layer in arcgis js api

texts:
I'm trying to zoom to layer automatically as soon as the map finishes loading using the guide here. To be clear, I'm only loading one layer and just want to zoom in on one of the layers provided by the service.

js script

dojo.require("dijit.layout.BorderContainer");
dojo.require("dijit.layout.ContentPane");
dojo.require("esri.map");
dojo.require("esri.renderer");
dojo.require("esri.graphic");
dojo.require("esri.dijit.Legend");
dojo.require("esri.utils");

dojo.require("dijit.layout.AccordionContainer");
dojo.require("esri.layers.FeatureLayer");

var map;
var basemap;
var trafficmap;
var basemapServiceSource;
var featuremapServiceSource;
var trafficServiceSource;
var renderer;
var loading;

function init() {

    // for now, these are referring to a single source
    basemapServiceSource = "http://co-gis-01/ArcGIS/rest/services/dev_RTIA/rtia/MapServer";

    //Add the topographic layer to the map. View the ArcGIS Online site for services http://arcgisonline/home/search.html?t=content&f=typekeywords:service
    basemap = new esri.layers.ArcGISDynamicMapServiceLayer(basemapServiceSource);
    basemap.setVisibleLayers([0, 1, 2, 4, 5, 6]);

    map = new esri.Map("map", {

            // remove logo
            logo:false
        });

    //loading image. id
    loading = dojo.byId("loadingImg");
    dojo.connect(map,"onUpdateStart",showLoading);
    dojo.connect(map,"onUpdateEnd",hideLoading);

    /*
        IDs of Layers for http://co-gis-01/ArcGIS/rest/services/dev_ExIS/exis/MapServer
            Site Name (0)
            Section ID (1)
            AADT some (2)
            AADT all (3)
            AdministrativeAreas_CongressionalDistricts (4)
            AdministrativeAreas_Provinces (5)
            AdministrativeAreas_Regions (6)
    */
    map.addLayer(basemap);

    // executes function on onLoad event
    dojo.connect(map, 'onLoad', function(theMap) {

        dojo.connect(dijit.byId('map'), 'resize', map,map.resize);
        zoomToLayer();
    });
}

function zoomToLayer() {

    var requestHandle = esri.request({
        url: "http://co-gis-01/ArcGIS/rest/services/dev_RTIA/rtia/MapServer/2",
        content: { f:"json" },
        callbackParamName: "callback"
    });

    requestHandle.then(requestSucceeded, requestFailed);
}

function requestSucceeded(response, io) {

    var extent = new esri.geometry.Extent(response.extent);
    dojo.byId("extent").innerHTML = dojo.toJson(extent.toJson());
    map.setExtent(extent);
}

function requestFailed()    {

    alert("requestFailed");
}

function showLoading()  {

    esri.show(loading);
    map.disableMapNavigation();
    map.hideZoomSlider();
}

function hideLoading(error) {

    esri.hide(loading);
    map.enableMapNavigation();
    map.showZoomSlider();
}

dojo.addOnLoad(init);

html code

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">

    <title>RTIA Map</title>
    <meta http-equiv="X-UA-Compatible" content="IE=EmulateIE8" >

    <link rel="stylesheet" type="text/css" href="gis.css" />

The problem is that the outputs of dojo.toJson(extent.toJson()) show the same extent, which is the full extent, no matter the layer. This particular layer that I would like to zoom to, http://co-gis-01/ArcGIS/rest/service...ia/MapServer/2, is just a filtered version of one of the "full-extent" layers. The filter is defined in the mapservice.

How do I correctly do this function?


texts: How do I drop the z value from the geometry and transfer it to a new field?

texts: Having a points shapefile with z values, how do I keep in the geometry only the x,y values? I also want to move the z values to a newly created double field.


texts: How to Zoom in QGIS Print Composer Without a Scroll Wheel

texts:
I am attempting to print a simple map from QGIS in portrait orientation.  When I use the "Add New Map" tool it brings the map in at an extent which matches my original data frame if I were to use landscape orientation.  See below:

However, when I set it to portrait, even if I do so before I add the map, it zooms out drastically and places my area of interest at the bottom end of the map.

Now, a bit of researching has revealed that I could zoom into the area I am interested in with my mouse's scroll wheel.  However, my mouse does not have a scroll wheel.  In addition, the Ctrl++ shortcut that works in the rest of QGIS does not work in Print Composer.  Any help would be appreciated.


texts: Compare two geometry tables using PostGIS

texts:
I have two geometries tables and I want to do a comparison between them two, so that I be able to select from one the geometries that don't exist  in the second. To do so I conceived the usage of  the following query :

create table touches_diammat_exists as select
touches_diammat.* from touches_diammat,pipe
where  st_astext(touches_diammat.fusion) != st_astext(pipe.geom)

Yet the results I'm getting have nothing to do with what I was expecting. It gave me a lot of geometries to the pace of a few. I don't know, I'm so confused and I even can't find  what I'm missing in my query.


texts: How do I get the points where a line crosses a polygon?

texts:
I have a line and need to get a list of coordinates where the line crosses a polygon. A line could cross multiple polygons. A line could also only cross one edge of polygon, if the line starts inside the polygon.

ST_CROSSES returns a bool indicating if the line crosses but I need the actual coordinates where the line intersects the polygon edge.

The basic SQL I have doesn't get me points:

SELECT trails.oid
FROM trails, areas
WHERE ST_Crosses(trails.geom, areas.geom);


texts: How to create a buffer considering the true surface area of the landscape?

texts: I would like to create a buffer around point. This is simple with ArcMap using the buffer tool. But is it possible to create the buffer considering also the true surface area of the landscape (consider also the slope distance using digital terrain model). In this case the buffer on slope will not be the circle but irregular ellipse. Are there any solution in ESRI or any packages in R?


texts:
I have a shapefile and I want to change the order of the features.
Is that possible?
I want to sort the features in the shapefile because I'm embedding this shapefile in a mobile app. The app will display a list of features, and although I could do the sorting at load time, I'd prefer to have them sorted.


texts: Digitise a circle around a point in QGIS 3?

texts:
In previous versions of QGIS I used the Intersect It plugin for drawing circles around a point with a specified radius.

However, I cannot see a way to digitise a circle around a point in QGIS 3.  I have tried using the new shape digitising tools, but these do not allow to centre the circle on a point (at least, not that I can see).

I have seen suggestions about buffering the points by a certain distance, but this seems quite a long-winded way of merely drawing a shape.

Any suggestions?


texts: Merging multiple bus stops with same name using QGIS?

texts:
My question is similar to Merging points (bus stops) in QGIS manually or programmatically?, but it's not answered sufficiently

I have OSM-transit data with bus stops. Some stops have two or more bays for different directions. In this example, the stop 'Steinerstraße' has two points:

I want to merge all points belonging to one stop to one point. I assume that they can be identified by an identical name in the 'name'-column of the attribute table.

Best case would be to create a new point in the middle of the points, however, I'd be satisfied with a simple deletion of the 2nd feature as well.

Can anyone give me a hint how to do it?

I've never used either python nor SQL in a GIS-environment and would prefer a GUI-based solution. However, I'm open for everything (some xp in python and SQL in other contexts).


texts: GeoServer - prioritize labels with certain attribute value

texts:
I have a point layer that I serve as WMS with GeoServer. The points have language and name properties and some points have the same geometry so they overlap. What I want to do is to prioritize rendering of labels for certain language using SLD stylesheet.

Example: There are two points with same geometry, one has language value 5, the other has language value 7. I want only the name of point with language value 7 to be rendered.


texts: Popup Not Displaying in ArcGIS API for JavaScript 3?

texts:
I would like to add a my popup work. I am confused which is which, infoWindow, popupTemplate and infoWindow. I have tried putting popupTemplate together but it does not load.
The basemap works alright but does not load the load the layer.

require([
"esri/map",
"esri/layers/FeatureLayer",
"esri/dijit/Popup",
"esri/dijit/PopupTemplate",
"esri/dijit/InfoTemplate",//new addition for infoTemplate
"dojo/domReady!"
 ],
 function(Map, FeatureLayer, Popup, PopupTemplate, InfoTemplate) {
//create map
var map = new Map("map", {
  basemap: "hybrid",
  center: [-0.345, 5.543],
  zoom: 19
});
//info template defination
var infoTemplate = new InfoTemplate("${OBJECTID}", "${*}");
var featureLayer = new FeatureLayer({
  url: "http://192.168.1.2/ndigis/rest/services/map/GSMA/FeatureServer/3",
  outFields: ["*"],
  infoTemplate: infoTemplate
});

map.addLayer(featureLayer);

});

based on the implementation here https://developers.arcgis.com/javascript/3/sandbox/sandbox.html?sample=renderer_class_breaks.


texts: Truncating tile cache outside bounding box using GeoWebCache

texts:
Is there a way to truncate all cached tiles outside of a specified bounding box, whilst leaving the cached tiles inside a bounding box untouched?
We have a nationwide system but it would be useful to remove all cached content outside of a bounding box whilst leaving the cached tiles inside the bounding box in situ. This would allow us to just reseed the cached content inside the bounding box when the mapping layers are updated.


texts: GDAL: Using gdallocationinfo on a rotated SAR image to get pixel value at Coordinates

texts:
Some information: I must not re-sample the image, the pixels are rotated, and the image opens in Arcmap just fine.
now, there are 15 GCPs on the raster (which make the pixels rotate)

Progress:
1.When I first try to use gdallocationinfo I get:
ERROR 6: No translation for an empty SRS to PROJ.4 format is known.

2.I use gdal_edit -ro -a_srs epsg:4326 "Image", this creates a .aux.xml file with the same name as my image, gdal now reads them both together, image still works on Arcmap.
then I use gdallocationinfo again and get:
ERROR 1: Cannot get geotransform

3.I use gcps2wld "Image" to get the geotransformation parameters, and insert it to the .aux.xml file via python GDAL dataset.SetGeoTransform

now Arcmap can't read the image.
and when I use gdalinfo, I get the following:

  Corner Coordinates:
  Upper Left  (   0.0000360,  -0.0000277) (  0d 0' 0.13"E,  0d 0' 0.10"S)
  Lower Left  (      -0.145,  887589.718) (  0d 8'43.59"W,Invalid angle)
  Upper Right (      -0.068,-1152505.040) (  0d 4' 3.91"W,Invalid angle)
  Lower Right (      -0.213, -264915.321) (  0d12'47.63"W,Invalid angle)
  Center      (      -0.107, -132457.661) (  0d 6'23.75"W,Invalid angle)

and if I Ignore that and use gdallocationinfo, coordinates that are supposed to be in the image give me errors for being outside the image borders.

I tried using gdal_edit -a_ullr to set the corner coordinates, but that just norths up the image (no rotation).

how do I set up the image right (without resampling) so I can get the right pixels in the right coordinates?


texts: Select the correct paths for The TIGER geocoder (linux)

texts:
I'm following this tutorial to setup the geocoder.

I'm stuck on the step 2.8.1.4

"edit the paths in the declare_sect column to those that fit Debbie's pg, unzip,shp2pgsql, psql, etc path locations."

The script is not finding any of the files.
E.G. for the PostgreSQL bin I have set:

export PGBIN=”/usr/lib/postgresql/9.3/bin”

Which is similar to the one in this answer:

And it seems to have the required files:

Nevertheless, I'm getting errors when I run nation_script_load.sh

nation_script_load.sh: 17: nation_script_load.sh: ”/usr/lib/postgresql/9.3/bin”/psql: not found

nation_script_load.sh: 18: nation_script_load.sh: ”/usr/lib/postgresql/9.3/bin”/psql: not found

Even the site to get the data seems to be falling:

ftp2.census.gov/geo/tiger/TIGER2015/STATE: No such file or directoryftp2.census.gov/geo/tiger/TIGER2015/STATE/.listing: No such file or directory


texts: Importing new layout to existing ArcMap file?

texts:
I am a new user of ArcMap and started with setting up data and I am now happy with the symbology styling I have set up. My data includes a mix of inhouse geodatabases, WMS feeds and shapefiles I create for specific one-off projects.

When I was being set up, another user set up a map layout which is not adequate for my purposes.  I have now located standard layout formats that I want to replace my current one with. Ideally I want to keep the general symbology styling the same, and change the layout to suit an individual project.

Is it possible to either replace/create a new map layout after one has already been started, or copy my data symbology into a new map based on the standard layouts?

I am using ArcGIS Desktop 10.7 on Windows 10.


texts: Select value from list in Python toolbox

texts:
within a Python toolbox in ArcGIS 10.2 I fill a parameter with a list. parameters[0].filter.list = myList It is essentialy a list of strings and will give the user a range of predefined years.
Later I check a file name if it holds information of a year. In case there is a year to be extracted from the name I want the value in the list defined for parameters[0] to be selected.

If I use SetParameterAsText() then I will overwrite my list of years, and cut my user short on the possibility to correct the year in case the extracted year is not the one the user is focusing on.

Since the year is anyway already part of the list... is there a way to say parameters[0].SelectValue() or parameters[0].SelectListItem() to get a visual confirmation for the user?

I hope it is understandable what I want to achieve.


texts: 'Select by Attribute' works outside of Graphic Modeler, but not within

texts:
I'm setting up a graphic modeler to handle three vector layers and one raster layer. In the end, I hope to vectorize the raster and merge the other three into one coverage. The raster file is a slope DEM, and I'm looking for anything steeper than a 25% grade.

I'm pulling my hair out because every time I use the 'select by expression' tool, the next step is to save the selection. Of course it tells me that nothing has been selected. I've tested my expression...

"Shape_Area">2000 AND "DN" = 1 (DN is the binary value attributed to slope < or > 25 degrees)

...outside of the graphic modeler and it works just fine. In debugging, I went back and saved the vectorized file with the $area field calculation included, and the 'select by expression' works just fine.

Thoughts?


texts: Making perpendicular lines through each polygon edge using PostGIS?

texts: I have polygons and lines I would like to make perpendicular lines through each polygon edge so I can find intersections with other polygons and lines.


texts: Displaying transparent overlapping polygons in ArcGIS for Desktop?

texts:
I'm trying to draw home range maps in ArcGIS 10.1 and am having difficulty displaying them as I would like.

The home ranges are in the form of polygons and many of them overlap. I'm trying to get a display setting that allows all polygons to be displayed at once, particularly where they overlap.

Essentially what I'm trying to do is display the polygons in the same way that the standard primary colour illustration appears, with three overlapping circles of different colours. The important part here is that where the polygons overlap they will combine to make a new colour, thus clearly demonstrating the overlapping area. Another way this is often done is to display each polygon with a unique crosshatching. For instance polygon 1 could have a left-facing diagonal hash, while polygon 2 could have a right-facing diagonal hash. Where they overlap, a cross-hatch pattern would be observed.

Transparency works in relation to other layers in the map, but unfortunately it doesn't seem to work in terms of each polygon within the layer. Same goes for defining the symbology level (essentially they all need to be on the same level).

I've seen this done with many other maps of animal home ranges or other overlapping zones, so I'm assuming there's got to be a way of doing this.


texts: Error when making a composite of four false color NDVIs

texts:
everything needed is here except importing Landsat 8 Surface Reflectance Tier 1.
I keep getting an error when I'm making a composite of 4 different false color NDVIs and I'm not sure what it means. I am somewhat of a beginner coding in GEE. Here is the error:

Here is my code:
var maskL8sr = function(image) {  // adapted from GEE docs
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3),
      cloudsBitMask = (1 << 5),
      snowBitMask = (1 << 4);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
               .and(qa.bitwiseAnd(cloudsBitMask).eq(0))
               .and(qa.bitwiseAnd(snowBitMask).eq(0));
  return image.updateMask(mask);
};
var addNdvi = function(image) { // adapted from GEE tutorial #6
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};
var favoriteBands = ['B2', 'B3', 'B4', 'B5', 'pixel_qa'];
dataset = dataset
    // .filterBounds(nj)  // this gets us paths we don't want
    .filterMetadata('WRS_PATH', 'equals', 14)
    .filterMetadata('WRS_ROW', 'equals', 32)
    .select(favoriteBands);
var
    winterFilter = ee.Filter.dayOfYear(356, 80),
    springFilter = ee.Filter.dayOfYear(81, 172),
    summerFilter = ee.Filter.dayOfYear(173, 264),
    fallFilter = ee.Filter.dayOfYear(265, 355),
    winter = dataset.filter(winterFilter).map(addNdvi).map(maskL8sr),
    spring = dataset.filter(springFilter).map(addNdvi).map(maskL8sr),
    summer = dataset.filter(summerFilter).map(addNdvi).map(maskL8sr),
    fall = dataset.filter(fallFilter).map(addNdvi).map(maskL8sr),
    diff = summer.select('NDVI').median()
        .subtract(winter.select('NDVI').median());

// set up some cartography stuff
// https://colorbrewer2.org/#type=sequential&scheme=Purples&n=5
var seqPalette = ['f2f0f7','cbc9e2','9e9ac8','756bb1','54278f'];
var rgbVis = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};
var colorizedVis = {
  bands: ['NDVI'],
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};
var diffVis = {
  min: -0.01,
  max: 0.46,
  palette: seqPalette,
};
states = ee.Image().byte().paint({
  featureCollection: states,
  width: 2,
});

// bake a layer cake
Map.centerObject(nj, 8);
Map.addLayer(winter, rgbVis, 'Winter RGB');
Map.addLayer(spring, rgbVis, 'Spring RGB');
Map.addLayer(summer, rgbVis, 'Summer RGB');
Map.addLayer(fall, rgbVis, 'Fall RGB');

Map.addLayer(winter, colorizedVis, 'Winter NDVI');
Map.addLayer(spring, colorizedVis, 'Spring NDVI');
Map.addLayer(summer, colorizedVis, 'Summer NDVI');
Map.addLayer(fall, colorizedVis, 'Fall NDVI');
Map.addLayer(diff, diffVis, 'Difference');
//Map.addLayer(states, {palette: ['DF0174']}, "States");

// debug
print(summer);
print(winter);
print(diff);
print(states);

// change detection RGB compositing

//Landsat 8
var vizParams1 = {
  bands: ['B5', 'B4', 'B3'],
  min: 0,
  max: 4000,
  gamma: [.9]
};

Map.setCenter(-74.0, 40.55, 8);
Map.addLayer(winter, vizParams1, 'T1 false color');
Map.addLayer(spring, vizParams1, 'T2 false color');
Map.addLayer(summer, vizParams1, 'T3 false color');
Map.addLayer(fall, vizParams1, 'T4 false color');
//concatenate 4 layers into a single image
var composite = ee.Image.cat([
  winter.select('B5'), spring.select('B5'), summer.select('B5'), fall.select('B5')]);
var vizParams2 = {
  bands: ['B5_2', 'B5_1', 'B5'],
  min: 0,
  max: 4000,
  gamma: [.9]
};
Map.addLayer(composite, vizParams2, 'T1-T2-T3-T4 NIR composite');


texts: How to display some contour lines and hide others?

texts:
I'm using ArcGIS 10.2 to show this contour lines information:

http://prntscr.com/4ncvc2

They're displayed with a 10 meter contour interval, but I want to display for example each 50 meters and hide those between lines. I'm using OS windows 7 and I'd like to get some heads up with this problem.


texts: Why is ArcPy / ArcGIS hanging when converting batches of large rasters to points?

texts:
There are 102 raster files in total that I want to convert to points. Each raster file is 2160x4320 pixels. When running the script in the code block below, the process hangs into the third or fourth raster file, regardless of if I take a random sample of the 102 rasters to start with. Even taking a subset of six of these rasters and converting them to points leads my machine to hang on the third raster file.

The first two rasters both take around 10 minutes to convert, and the third raster file was processing for 10 hours before I stopped it. While the third raster file was being processed, I noticed that ArcMap's associated RuntimeLocalServer.exe private memory allocation for RAM grew from ~900mb after completing the first two and starting on the third to ~2600mb once I stopped the process.

Regarding the machine I'm using - I'm using ArcGIS 10.5 on a computer with a i7-6700K processor, 32GB ram, and a terabyte SSD. The computer is good working order in terms of software and hardware.

When run on a desktop a few years older with only 16GB of RAM, the same script below has crashed ArcGIS and left the routine incomplete after about three files are processed. It seems like ArcGIS isn't clearing VRAM or RAM used for already-converted raster files, but I have no idea if this is actually the problem.

Is there some ArcPy memory management command that I can include in my code below? Why is this process consistently failing to finish?

# Convert all rasters in a directory to points.

import arcpy
import os

conversionList = []

for root, dirs, files in os.walk(r"D:/data/crops/earthstat/HarvestedAreaYieldMajorCrops_Geotiff/"):
    for f in files:
        if f.endswith(".tif"):
            fpath = root + "/"+  f
            conversionList.append(fpath)

outputPath = r"D:/data/crops/earthstat/HarvestedAreaYieldMajorCrops_Points/"

for item in conversionList:
    infile = str(item)
    outfile = outputPath + item.split("/")[-1].split(".")[0] + ".shp"
    arcpy.RasterToPoint_conversion(infile, outfile, "Value")


texts: How to determine why feature is invalid according to Validate features?

texts:
I have a feature class with a lot of domains and a lot of different properties e.g. null values, subtypes, and so on.

I am using Validate features (Editor tool) and get the selection of the invalid feature(s), but I cannot figure out why some features are invalid.

Is there a way to get a report or something?


texts: Keep the popover in the same position when zooming in the map with Bootstrap 4 and OpenLayers 5

texts:
I tried the popover example from the OpenLayers web site: https://openlayers.org/en/latest/examples/overlay.html

It works very well but when I change the version of Bootstrap to version 4 and jQuery to version 3.3.1 the popover change the location when I zoom on the map.

I want to keep the popover in the position that I clicked even when I zoom.


texts:
I am researching about shapefiles. I am not able to understand the usage of polyline Z parameters.

I have searched the internet but I didn't come across any example where this parameter is being used.

Please help me to understand this parameter better?


texts: I have a list of geographic locators and a boundary file. In QGIS, is there any way merging them together to create a lookup saying this boundary contains these locators? Not the number the actual locator i.e. postcode, easting/northing.


texts: Configurable Map Viewer (CMV) Measurement Widget Error

texts:
I have several cmv widgets enabled, some I have built myself. All work great, but when I enable the Measurement widget I am getting this error upon loading the app:

  dojo/parser::parse() error TypeError: Cannot read property 'toString'
  of undefined

  Uncaught Error: uniqName_4: parser returned unfilled promise (probably
  waiting for module auto-load), unsupported by _WidgetsInTemplateMixin.
  Must pre-load all supporting widgets before instantiation.

Here is my configuration:

  measure: {
    include: true,
    id: 'measurement',
    type: 'titlePane',
    canFloat: true,
    path: 'gis/dijit/Measurement',
    title: ' &nbsp;&nbsp; Measure',
    open: false,
    position: 4,
    options: {
      map: true,
      mapClickMode: true,
      defaultAreaUnit: units.ACRES,
      defaultLengthUnit: units.FEET
    }
  }

Any ideas what the problem might be?


texts:
Here is a Google Earth screenshot of “Edwards AFB” in California:

-The green circle has the compass in it...
Here’s a zoomed in screenshot of the compass:

Notice how my Google Earth compass (in green again) is pointing directly North but the compass in the picture isn’t very close to pointing North compared to the other one
If I had to guess which compass is off, I would say it’s Google’s, but why wouldn’t the compasses be closer than this?


texts: Dissolving polygon features with duplicate set of attributes

texts:
I have one shapefile with thousands of polygon features. Many of them have the exact sets of attributes. So I want to make those features into one so that I can reduce the number of features from a few thousands to, say, ten to twenty.

I found the Polygon dissolve (by attributes) tool quite similar to what i want but I can only input 3 attributes in the GUI. But I want to input all attributes. Is it possible to do it in PyQGIS?


texts: How can I change layerswitcher properties(radio button) ? So I want to change just like drop down list ?


texts: What's the utility of osm_meta column in osm2po (Version: 4.6.9) generated data?


texts: QGIS: Filtering aggregate returns different results with same value

texts:
I have a joined layer that was created from a one-to-many join between a vector layer with three attributes (lat, long, siteid) and a data layer that contains many columns with many entries possessing identical 'siteid' fields. I need to find the minimum and maximum range of a field called 'sampleno' for each group of entries containing the same 'siteid'.

Here's what I tried:
maximum(sampleno,'', siteid=attribute($currentfeature,'siteid'))
|| '\n' ||
maximum(sampleno,'',siteid='1750')

I also tried:
 maximum(sampleno,'', siteid="siteid"))
 || '\n' ||
 maximum(sampleno,'',siteid='1750')

I don't understand why entering a test value of '1750' as the filter returns the result I expect but a value of "siteid" or attribute($currentfeature,'siteid') doesn't.


texts:
Is it possible to rename the favourites in QGis Browser ? My folder paths are pretty long and I'd like to see directly my folder name.

I'm sorry if it's already been answered, I've looked online but found nothing on this simple question.


texts:
I've built a model that takes the name of the Excel worksheet containing xy data, and uses that name as the name of the output shapefile. ( Worksheet "massachusetts_sightings" in an xlsx becomes "massachusetts_sightings$.shp"). I do this using Parse Path. Notice that a "$" is added.

If I then feed this shapefile, massachusetts_sightings$.shp into a clip tool (say, to clip out points in the data that are not in my study area), I get an error for an invalid character in the clip tool. My suspicion is the $. I know I could remove this with Calculate Value...but I can't figure out the code block.

My parse path output is called Value:

Can someone help troubleshoot? I haven't used Python in years-

# Import arcpy module
import arcpy

# Load required toolboxes
arcpy.ImportToolbox("Model Functions")

# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "X:\\Tools\\XY\\point_upload\\Scratch"
arcpy.env.workspace = "X:\\Tools\\XY\\point_upload\\Data"

# Script arguments
XY_Table = arcpy.GetParameterAsText(0)

Define_Coordinate_System = arcpy.GetParameterAsText(1)

Projection = arcpy.GetParameterAsText(2)

Geographic_Transformation = arcpy.GetParameterAsText(3)

Study_Area = arcpy.GetParameterAsText(4)

Projected_and_Clipped_Data = arcpy.GetParameterAsText(5)
if Projected_and_Clipped_Data == '#' or not Projected_and_Clipped_Data:
    Projected_and_Clipped_Data = "X:\\Tools\\XY\\point_upload\\Data\\%fixed%_Clip.shp" # provide a default value if unspecified

# Local variables:
xy_event = XY_Table
v_defined = xy_event
defined_data = v_defined
v_projected = defined_data
Value = XY_Table

# Process: Parse Path
arcpy.ParsePath_mb(XY_Table, "NAME")

# Process: Make XY Event Layer
tempEnvironment0 = gp.scratchWorkspace
arcpy.env.scratchWorkspace = "X:\\Tools\\XY\\point_upload\\Scratch"
tempEnvironment1 = gp.workspace
arcpy.env.workspace = "X:\\Tools\\XY\\point_upload\\Data"
arcpy.MakeXYEventLayer_management(XY_Table, "", "", xy_event, "", "")
arcpy.env.scratchWorkspace = tempEnvironment0
arcpy.env.workspace = tempEnvironment1

# Process: Copy Features
tempEnvironment0 = gp.scratchWorkspace
arcpy.env.scratchWorkspace = "X:\\Tools\\XY\\point_upload\\Scratch"
arcpy.CopyFeatures_management(xy_event, v_defined, "", "0", "0", "0")
arcpy.env.scratchWorkspace = tempEnvironment0

# Process: Define Projection
arcpy.DefineProjection_management(v_defined, Define_Coordinate_System)

# Process: Project
tempEnvironment0 = gp.scratchWorkspace
arcpy.env.scratchWorkspace = "X:\\Tools\\XY\\point_upload\\Scratch"
tempEnvironment1 = gp.workspace
arcpy.env.workspace = "X:\\Tools\\XY\\point_upload\\Data"
arcpy.Project_management(defined_data, v_projected, Projection, Geographic_Transformation, "")
arcpy.env.scratchWorkspace = tempEnvironment0
arcpy.env.workspace = tempEnvironment1

# Process: Clip
arcpy.Clip_analysis(v_projected, Study_Area, Projected_and_Clipped_Data, "")

# Process: Calculate Value
arcpy.CalculateValue_management("\"%Value%\".strip(\"$\")", "\\n", "String")


texts:
I am new to GeoServer. I have a PostGIS database and GeoServer setup. How can I query the database using the GeoServer/ Boundless SDK?

I have OpenGeo Suite 4.6 and I am trying to query the database in a custom WPS process. I don't want to directly connect to the database but use GeoServer SDK somehow to query the db.


texts: How do I install PostgreSQL + PostGIS + MapServer + Python-MapScript on Ubuntu 12.04 LTS w/o compiling from source?

texts:
This question is similar to How to have PostGIS 2.1 and the Python osgeo library installed together without compiling from source?

but I couldn't figure out how to ask a follow-up question.

I am trying to get PostgreSQL + PostGIS + MapServer + Python-MapScript working on Ubuntu 12.04 LTS, for a Geo-Django site. In the past, this has worked wonderfully (was still working in April of 2014), however something has changed, and I'm now getting an error regarding gdal compatibility:

The following packages have unmet dependencies: postgresql-9.1-postgis : Depends: libgdal1 (>= 1.9.0) but it is not going to be installed E: Unable to correct problems, you have held broken packages.

After a fair bit of searching, I came across Brian (cokrzys)'s answer to the post above. I followed his instructions on a clean install of Ubuntu 12.04, and was subsequently able to create the PostGIS templates in PostgreSQL 9.3. However, when I then try to install python-mapscript, I get the following message:

  The following packages will be REMOVED:
   libgdal1 postgis postgresql-9.3-postgis-2.1

So, the problem is that it does not appear that I can have both the postgis libraries and the python-mapscript libraries installed at the same time.

Honestly, I don't care what version of PostgreSQL / PostGIS / MapServer I use, so long as I can get it running on Ubuntu 12.04.

Has anyone had any luck with this?

@Mark helped me to finally got this working as follows:

#!/usr/bin/env bash
###################################################################
# Install GDAL, MapServer, etc. first
###################################################################
echo "Y" | sudo apt-get install python-software-properties
echo "Y" | sudo add-apt-repository ppa:ubuntugis/ppa
sudo apt-get update
echo "Y" | sudo apt-get install mapserver-bin
echo "Y" | sudo apt-get install gdal-bin
echo "Y" | sudo apt-get install cgi-mapserver
echo "Y" | sudo apt-get install python-gdal
echo "Y" | sudo apt-get install python-mapscript
echo "Y" | sudo apt-get install python-gdal

###################################################################
# Then Install PostgreSQL9.1, PostGIS 9.1
###################################################################
echo "Y" | sudo apt-get install postgresql-9.1
echo "Y" | sudo apt-get install postgresql-client-9.1
echo "Y" | sudo apt-get install postgresql-server-dev-9.1
echo "Y" | sudo apt-get install postgresql-plperl-9.1
echo "Y" | sudo apt-get install postgresql-9.1-postgis-2.0


texts: Using CopyFeatures_management without copying source GDB domains to 'in_memory' workspace?

texts:
I know it's probably the opposite of what would generally be desirable, but I am trying to temporarily copy a Feature Class from a File GDB with Domains into a in_memory intermediate FC without maintaining the source Domains.

I cannot prove this because I don't know how (where) to look for the in-memory FC I produce, but when I do

temp_folder = r'in_memory'
arcpy.arcpy.CopyFeatures_management(fc, os.path.join(temp_folder, 'temp'))

,I think I not only copy the entire source FC into the in_memory workspace, but also the source FileGDB's Domains.

Can someone confirm this behaviour? Is there a way to avoid this using the same tool (CopyFeatures_management) or (only if it's not possible) a similar tool?

I would like to keep using in_memory for my intermediate file, rather than create physical intermediate files because:

I would have to delete them
I would have to properly check for errors within my script, handle them in order to delete intermediate files before they reach the 'delete' part of the script
in_memory is faster because it uses RAM (is this true?)

I tried deleting the domains, but the code below seems to make ArcGIS crash all the time (imagine domains comes from something like arcpy.ListDomains)...

for d in domains:
    arcpy.DeleteDomain_management(temp_folder, d)


texts: SQL MERGE on Remote server - Geometry Data Type Objects exposing columns with CLR types are not allowed in distributed queries

texts:
I'm trying to merge tables. Both contain a data type of geometry. Apparently I need a passthrough query to be able to write geometry/geography type to remote server.
MERGE
  [SERVER].[DB].[dbo].[TABLE] AS TR
USING
  #Temptable AS SR
ON
  (TR.ASSET_ID = SR.ASSET_ID)

WHEN MATCHED AND TR.SHAPE <> SR.SHAPE OR TR.VERS <> SR.VERS
    THEN UPDATE SET
        TR.SHAPE = SR.SHAPE,
        TR.VERS = SR.VERS

WHEN NOT MATCHED BY TARGET
    THEN INSERT (ASSET_ID, VERS, SHAPE)
        VALUES (SR.ASSET_ID, SR.VERS, SR.SHAPE);

Error in SQL I'm getting:

Objects exposing columns with CLR types are not allowed in distributed queries. Please use a pass-through query to access remote object


texts: PostGIS cannot create extention:ERROR: cache lookup failed for function 1

texts:
When trying to create PostGIS extension I am getting this error:

  ERROR: cache lookup failed for function 1

I am running PostgreSQL 11 on Ubuntu 16.04.
I have installed PostGIS and PostGIS scripts. I managed to create extension on one database but something caused this error so that I can't create it anymore.


texts: Problem in splitting and removing Polyline

texts:
I am stuck at a problem for weeks, Please Help!!!

I am making program in C#
Here is what I need to do --- I have two polylines that are intersecting with each other at some points. I want to delete those portion of one polyline which are going above the other polyline.
I am attaching an image explaining the problem I am going through. In this image, I have two Polylines Line1 and Line2. So according to that image, I want to remove the dashed portions of Line2.

Please Help!!!
Thanks a ton in advance

Here is what I have done till now....
First I am finding out all the points of intersection of the two lines and saving them in a layer....
Then I am Splitting these lines at these points using this code...

while (pPointF != null)
{
    IPoint pPoint = (IPoint)pPointF.Shape;
    ISpatialFilter pSF = new SpatialFilterClass();

    pSF.Geometry = pPoint;
    pSF.GeometryField = "Shape";
    pSF.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;

    IFeatureCursor pLineCursor = pFeatureclass1.Search(pSF, true);
    IFeature pLineF = pLineCursor.NextFeature();

    while (pLineF != null)
    {
        IPolycurve pPolycurve = (IPolycurve)pLineF.Shape;
        IPoint pToPoint = pPolycurve.ToPoint;
        IPoint pFromPoint = pPolycurve.FromPoint;

        if (pFromPoint.X != pPoint.X && pFromPoint.Y != pPoint.Y && pToPoint.X != pPoint.X && pToPoint.Y != pPoint.Y)
        {
            IFeatureEdit pFeatureEdit = (IFeatureEdit)pLineF;
            pFeatureEdit.Split(pPoint);
        }
        pLineF = pLineCursor.NextFeature();
    }
}

The problem with this code that it is giving different results than expected.
As I can see from attribute table, it is generating more than expected parts of line i.e. if there are 5 intersection points then there should be 6 parts of a line, nut it is creating 10-12 parts. Also the number of parts created for both lines are different, which logically should be same.
Due to these problems, when I am trying to delete unwanted parts of line, the arcGis crashes.....
Please help with this problem.
.
.
.
Now I am using a Geoprocessor to do splitting, following is the code:

        SplitLineAtPoint pSplitter = new SplitLineAtPoint();
        pSplitter.in_features = pPolylayer1;
        pSplitter.point_features = pPointsFeatureLayer;
        string outfile = txtDirName.Text + "\\" + "splitted";
        pSplitter.out_feature_class = @outfile;

        Geoprocessor GeoPro = new Geoprocessor();
        GeoPro.Execute(pSplitter, null);

But on the last line, where Geoprocessor is being executed, the ArcGIS crashes...
Please tell me where is the problem...


texts: GeoServer WMS GetMap request from multiple layers

texts:
I'm trying to form a GetMap request for my geoserver that returns data (.png image) from two layers (buildings and roads).

Request looks like this:

https://mysite.net/geoserver/wms?request=GetMap&service=WMS&version=1.1.1&layers=buildings,roads&styles=buildings,roads&srs=EPSG%3A3765&bbox=397034,4909706,397213,4909832&&width=600&height=600&format=image%2Fpng

And this is the response:

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE ServiceExceptionReport SYSTEM
"https://mysite.net/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd">
<ServiceExceptionReport version="1.1.1" >
    <ServiceException code="LayerNotDefined" locator="layers">
      Could not find layer  roads
</ServiceException>
</ServiceExceptionReport>

If I want data from only one layer, it works. Like example below (returns image with buildings only):

https://mysite.net/geoserver/wms?request=GetMap&service=WMS&version=1.1.1&layers=buildings&styles=buildings&srs=EPSG%3A3765&bbox=397034,4909706,397213,4909832&&width=600&height=600&format=image%2Fpng

How can I request an image that contains both layers?


texts: is it possible to convert 9.3 lyr to 10 lyr

texts:
I have a customer who sent me a lyr file in 9.3 and I am using 10,  is there a way to import it in.  It seems to hang up whenever I attempt to open it!

Thanks


texts: QGIS save postgis table error : no result buffer

texts:
I've a polygon layer on PostGIS, and when I try to add featureS using QGIS it gets error when I save the edit.

Here is the error notification "no result buffer".

I found a similar case here : PostGIS error while adding feature: no result buffer. But it doesn't give me the answer to solve it. Any ideas ?


texts: How do I insert a date string into the database as a date?

texts:
I'm using an InsertCursor to insert rows into a feature class. One of the fields in the source data may contain a NULL data value. The destination feature class requires a date - it is not nullable. In the dictionary of values I pass to the InsertCursor, I insert a string value of "01/01/1970 00:00:00". When the cursor tries to insert the row, I get an error that says the column cannot be null.

How do I get the InsertCursor to insert the date string as a date?


texts: Changing the colour of a GeoJSON polygon if user location is within it in Leaflet

texts:
I am wondering if there is a way to change the colour of a polygon got from a geojson file only if the users/browser location is within that polygon if so is there any example or working code that cpild be shared
No code But some pseudocode to give the general idea:

Function locationWithinPolygon (e) {  if (browserlocation within
GeoJSONPolygonBoundaries){   pan to feature  Change colour of
polygon to white   popUpMessage "Your are within" + polygonName }
 else {  function OnLocationError{ alert(e.message); }
map.on("location error" + OnLocationError)   }


texts:
Short version

Turn restrictions don't seem to be taken into account in GraphHopper's routing, although their GitHub issues say it is implemented. Does anybody know the state of this? Is it implemented, but not merged, not released, not rolled out or simply not turned on? Or is there a problem with the data in the example below? What's the state of turn restrictions in other open OpenStreetmap-based routing solutions?

Detailed version

Situation

GraphHopper Maps is my go-to solution for routing. It seems the best routing service with OpenStreetmap data, which in turn is the most complete database for cycling ways and footpaths.

Since I'm also doing some OpenStreetmap mapping and wanted to advance to more sophisticated features like turn restrictions, I tried to learn from other people's examples. They seem to be recognized by JOSM, but not taken into account by GraphHopper.

Example

Here an example of medium complexity with a crossing that contains a right-turn restriction and a U-turn restriction. OpenStreetmap seems to have the proper turn restrictions and JOSM seems to understand them (judging from the traffic signs shown in the graphical representation).

Note that I selected the "from" via in both cases, so it gets selected in the map visualization and you can confirm that there is no confusion between the "from" and "to" entry.

But when I try it on GraphHopper it turns out that neither of the two restrictions are taken into account. Both the right-turn and the subsequent U-turn in the route that GraphHopper proposes are actually illegal:

GraphHopper implementation

So I was wondering, if the problem is in the OpenStreetmap data or in GraphHopper.

From the GitHub issue #2 "Support for turn restrictions and costs" I understand that the feature of turn restrictions was actually implemented a while ago. But there is a related GitHub issue #55 "Initial implementation of turn restriction support by khuebner" that talks about some merging problems with too complex code in a single commit. It also references the merge being split in several parts, the pull requests #133, #134, #135 and #158. They are all marked as "closed" or "merged".

Questions

So does anybody know what's the state of this? Is it implemented, but not merged, not released, not rolled out or simply not turned on? Or is there a problem with the data in this example? What's the state of turn restrictions in other open OpenStreetmap-based routing solutions?


texts: OpenLayers6 TileArcGISRest LAYERS param doesn't allow to modify and transparency of the features

texts:
I am trying to modify the transparency of the features ans I don't have a clue about how to do it.

The site https://openlayers.org/en/latest/apidoc/module-ol_source_TileArcGISRest.html refers to https://developers.arcgis.com/rest/services-reference/export-map.htm and I don't know is it possible, I tried this:

params: {
  'layers':"show:3",
  'drawingInfo': "transparency: 50",
},

but transparency could be an attribute of the layer.

Any idea?


texts: I've been trying without success to find free textures to use in TileMill. I'm looking for fairly simple, abstract textures, like a dotted pattern (to indicate national parks). Are there any well known resources for this kind of thing?


texts: QGIS 3.2 - number field display in attribute table in american format

texts:
Since upgrading to QGIS 3.2.0 the int and numeric formatted fields are displaying in the attribute table with a comma added for numbers of more than 3 digits (e.g. 12,345 instead of 12345).

It looks like I am in an American version and I assume it is just a setting but I have searched through the options/settings/properties and no success so far.

I have also reinstalled QGIS 3.2 to make sure I haven't ticked an American installation option but no change.

How can I 'switch off' that comma??

I previously run QGIS 2.18.13 and the number display in the attribute table was without the comma.


texts: ArcGIS for Desktop adds fields to *.tab table?

texts:
I have a tab-delimited table. The file is a .txt file with the .tab extension. The rows are:

Dataset Feature_Class_or_Table  Field   GLOBALID    New_Value   Bad_Character
GasDataset  CasingAnno  TEXTSTRING  {5A5B2C4B-F395-4AD2-BEB8-D1AB3D61DDC2}  3/4 -4  new line, carriage return
GasDataset  CasingAnno  TEXTSTRING  {832D9957-CEA2-4DDE-801A-4A12F286681C}  2 PL-34 new line, carriage return
GasDataset  CasingAnno  TEXTSTRING  {66441560-5EF1-4572-B13B-F7C4970CD2C1}  2 PL-30 new line, carriage return
GasDataset  CasingAnno  TEXTSTRING  {EF2B66AF-E957-47AA-A23F-1334E085C5A4}  1 1/2 PL    new line, carriage return

However, when ArcGIS 10.1 for Desktop reads the table, it adds a couple new fields; New_Value_X and New_Value_Y.

Is there a way to prevent the reading/creation of these two fields?


texts: I've got a .csv file that contains a long list of latitude and longitude points. Is there an efficient way to batch convert the list into a projected UTM coordinate system?


texts: TileCache SVN trunk on IIS 7 - An error occurred: No section: 'cache'

texts:
I was having tile generation problems with TileCache 2.11 and decided to update to the trunk.  I repeated the instructions necessary to get TileCache working in IIS but now I'm seeing this error:

An error occurred: No section: 'cache'
  File "C:\inetpub\wwwroot\tilecache\TileCache\Configs\File.py", line 147, in read
    self.cache = self._loadFromSection(config, "cache", TileCache.Cache)
  File "C:\inetpub\wwwroot\tilecache\TileCache\Config.py", line 74, in _loadFromSection
    type  = config.get(section, "type")
  File "C:\Python25\lib\ConfigParser.py", line 511, in get
    raise NoSectionError(section)

I've made the necessary change to tilecache.py:

from TileCache import Service, cgiHandler, cfgfiles

if __name__ == '__main__':
    svc = Service.load("c:\\inetpub\\wwwroot\\tilecache\\tilecache.cfg")
    cgiHandler(svc)

And I've confirmed my configuration file exists and does contain a [cache] section:

[cache]
type=Disk
#base=c:\inetpub\wwwroot\tilecache\Cache
base=g:\tilecache\Cache

Using SysInternals ProcessMonitor I've determined that python.exe never attempts to access tilecache.cfg.

Renaming tilecache.cfg results in the same error.  Why isn't it loading the file?  Is the path incorrectly formed?


texts: "User memory limit exceeded" during Random Forest classification in GEE

texts:
I succeeded in an initial attempt at using Random Forest classification on my input image (an exported mosaic covering ~ 2,797 km2). I used 30 trees, classified based on 4 bands, and passed a training set that included 251 polygons that covered ~3 km2. I wanted to refine the outputs, so I added training polygons in an attempt to have the training set cover 0.25% of the total study area (65 polygons per class @ 1.4 km2 each ~= 7 km2).
After boosting my training set and re-importing, I now get a "User memory limit exceeded" error. I have tried adding a tileScale property to my image.sampleRegions() parameters, but this does not solve the issue:
var imageName = ee.Image(".../201808image");

//Set visualization parameters for PlanetScope tiles
var vizParams = {
  bands: ['b3', 'b2', 'b1'],
  min: 250,
  max: 1452,
  gamma: 1.0
};

//Set up supervised classification - use all PlanetScope bands.
var bands = ['b1', 'b2', 'b3', 'b4'];
var p = [
  '369b47', //forest
  'b76031', //low veg
  'cc0013', //built up
  'cdb33b', //tilled land
  'aec3d4' //water
];

//Load training dataset
var trainSet = ee.FeatureCollection('.../TrainingSet');

//Overlay polygons and begin training
var training = imageName.sampleRegions({
  collection: trainSet,
  properties: ['Class'],
  scale: 3,
  tileScale: 16
});

var classifier = ee.Classifier.smileRandomForest(30).train({
  features: training,
  classProperty: 'Class',
  inputProperties: bands
  });
print('Random forest, explained', classifier.explain());

// Classify the image with the same bands used for training.
var classified = imageName.classify(classifier);
print('Image classified!');

//Add image created from PlanetScope tiles to map
Map.addLayer(imageName, vizParams, 'PlanetScope imagery 201808');
print('Image added.');
Map.addLayer(trainSet, {color: 'aqua'}, 'Training Polygons');
print('Training polygons dataset added.');
Map.addLayer(classified, {palette: p, min: 0, max: 4}, 'land cover');
print('Classified layer added.');

How do I get around the user memory limit that I am exceeding?
I get the impression people use GEE to process much larger image datasets. I am unclear if there is something with my training polygon Feature class or in the creation of SampleRegions that I should be changing.


texts: Finding closest Intersection between line strings (PostGIS)

texts:
I want to find the distance between one point (excavation sites) and a line (rivers) for the cardinal directions.
What I thought would be the easiest way is to draw lines between the excavation site and a given maximum distance every 45° and find the intersection between the line and the river which isn't working properly yet (and maybe isn't the best idea to solve that?).

I have a table excavation_sites and a table dist_ex_sites_river (in which the lines and the distances should be stored). The dist_ex_sites_river table also has a foreign key for the excavation sites. I managed to create the lines already and write them in the db with the according step_nr (see figure - the red circles are the intersection points I want to find for measuring the distance) and the excavation site. Now I'd like to get the intersections to measure the distance which always returns None, maybe you guys can help me?
This is the code I'm using right now:

SELECT ST_Intersection(dist_ex_sites_river.line_geometry,env_white_water.geom) FROM public.env_white_water, public.dist_ex_sites_river WHERE dist_ex_sites_river.step_nr_id_pk = 1 AND env_river.river_id_pk = 616;

I'm working with python and the psycopg2 package. I iterate through every id (in both tables) and want to append the distances to a list and choosing the minimum value in the end to insert that into the db in the end.
Probably this is absolutely stupid and there is a much better way to do that.
Thanks a lot for your help.


texts: How to connect the selectionChanged signal with pyqt?

texts:
I am developing a plugin for QGIS in Python, I need to connect the sselectionChanged signal emitted when a feature of the layer is selected, I could not find any examples on internet, here is what I have done so far:

QObject.connect(self.iface.mapCanvas(),SIGNAL("selectionChanged(QgsMapLayer)"), self.test)

In the logs of QGIS I get this warning:

Warning: Object::connect: (receiver name: 'MainWindow')

can you help me?


texts: Extract USGS satellite image inside a given bounding box

texts:
I'm using USGS explorer to extract satellite image. What I want is that given a bounding box, I want to extract satellite image inside the box (see the image attached). I don't want the blur satellite images, instead, I want to get the clear image with same quality as in the USGS explorer tool.
Is there a programming way (using Python or any programming language) to obtain these images with USGS explorer or any other data source?


texts: Problem when I try to a time series chart Error generating chart: No features contain non-null values of "system:time_start"

texts:
I've calculated NDVI from Sentinel 2 images but when I try to get the chart of the series this message is generated - Error generating chart: No features contain non-null values of "system:time_start".
 The code is the following:

var s2 = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2019-01-01', '2019-12-31')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .filterBounds(studyArea);

function getNDVI (image) {
    var ndvi = image.normalizedDifference(['B8','B4']).rename("ndvi");
    image = image.addBands(ndvi);
    return image;
}
s2 = s2.median();
s2 = getNDVI(s2).clip(studyArea);
var s2NDVI = s2.select('ndvi');
print(s2NDVI);

var NDVIseries = ui.Chart.image.seriesByRegion(
    s2NDVI, studyArea, ee.Reducer.mean(), 'ndvi', 200,  'system:time_start', 'label')
        .setChartType('ScatterChart')
        .setOptions({
          title: 'NDVI 2016',
          vAxis: {title: '-1+1'},
          lineWidth: 1,
          pointSize: 4,
          series: {
            0: {color: 'FF0000'},
}});

// Display.
print(NDVIseries);

Can somebody help me to fix the code? Thanks a lot.

Paolo


texts: Nominatim: osm2pgsql path specification problem

texts:
The Problem:

I've installed Nominatim on CentOS 7, but I'm running into an import error when I try to load the data. Here is the error:

Import
Please download and build osm2pgsql.
If it is already installed, check the path in your local settings  (settings/local.php) file.
ERROR: osm2pgsql not found in '/opt/Nominatim/osm2pgsql/osm2pgsql'
osm2pgsql not found in '/opt/Nominatim/osm2pgsql/osm2pgsql'

Osm2pgsql is already installed and here are the following paths after I run some find commands:

/usr/local/bin/osm2pgsql
/usr/local/share/osm2pgsql

I've tried specifying these paths in the local.php file like this:

@define('CONST_Osm2pgsql', '/usr/local/bin/osm2pgsql');
@define('CONST_Osm2pgsql', '/usr/local/share/osm2pgsql');

Unfortunately, I still get the same problem, so I must not be specifying the paths right. Does anyone know the syntax of local.php? There aren't any guides I could find online beyond basic postgresql and postgis configurations.


texts: Can satellite data have different locations of nodata values in different channels

texts:
I have some Sentinel satellite imagery.

And I have noticed that it looks like there are different locations of nodata values in each channel's mask. (I read in the tif file with 3 channels in via the rasterio python library)

This seems VERY strange to me. Shouldn't the nodata masks be the exact same across channels?

This file is a Sentinel 2 Infrared Ratio band composite image (3 channels) (the yellow dots are the nodata pixels):

I don't have a background in satellite imagery / GIS though, so perhaps this actually is normal? Could someone explain to me how this could be possible.

Or for satellite data, should I only consider a pixel as nodata if ALL 3 bands have the nodata value? Because at those locations where one channels show nodata and others don't, the values at all 3 channels are for example [3, 100, 255] where 255 is nodata.


texts:
A common requirement in GIS is to apply a processing tool to a number of files or apply a process for a number of features in one file to another file.

Much of these operations are embarrassingly parallel in that the results of the calculations in no way influence any other operation in the loop. Not only that but often the input files are each distinct.

A classic case in point is the tiling out of shape files against files containing polygons to clip them against.

Here is a (tested) classical procedural method to achieve this in a python script for QGIS. (fyi the output of temporary memory files to real files more than halved the time to process my test files)

import processing
import os
input_file="/path/to/input_file.shp"
clip_polygons_file="/path/to/polygon_file.shp"
output_folder="/tmp/test/"
input_layer = QgsVectorLayer(input_file, "input file", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(input_layer)
tile_layer  = QgsVectorLayer(clip_polygons_file, "clip_polys", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(tile_layer)
tile_layer_dp=input_layer.dataProvider()
EPSG_code=int(tile_layer_dp.crs().authid().split(":")[1])
tile_no=0
clipping_polygons = tile_layer.getFeatures()
for clipping_polygon in clipping_polygons:
    print "Tile no: "+str(tile_no)
    tile_no+=1
    geom = clipping_polygon.geometry()
    clip_layer=QgsVectorLayer("Polygon?crs=epsg:"+str(EPSG_code)+\
    "&field=id:integer&index=yes","clip_polygon", "memory")
    clip_layer_dp = clip_layer.dataProvider()
    clip_layer.startEditing()
    clip_layer_feature = QgsFeature()
    clip_layer_feature.setGeometry(geom)
    (res, outFeats) = clip_layer_dp.addFeatures([clip_layer_feature])
    clip_layer.commitChanges()
    clip_file = os.path.join(output_folder,"tile_"+str(tile_no)+".shp")
    write_error = QgsVectorFileWriter.writeAsVectorFormat(clip_layer, \
    clip_file, "system", \
    QgsCoordinateReferenceSystem(EPSG_code), "ESRI Shapefile")
    QgsMapLayerRegistry.instance().addMapLayer(clip_layer)
    output_file = os.path.join(output_folder,str(tile_no)+".shp")
    processing.runalg("qgis:clip", input_file, clip_file, output_file)
    QgsMapLayerRegistry.instance().removeMapLayer(clip_layer.id())

This would be fine except that my input file is 2GB and the polygon clipping file contains 400+ polygons. The resulting process takes over a week on my quad core machine. All the while three cores are just idling.

The solution I have in my head is to export the process out to script files and run them asynchronously using gnu parallel for example. However it seems a shame to have to drop out of QGIS into an OS specific solution rather than use something native to QGIS python. So my question is:

Can I parallelise embarrassingly parallel geographic operations natively inside python QGIS?

If not, then perhaps someone already has the code to send this sort of work off to asynchronous shell scripts?


texts: Changing kml marker in esri javascript api

texts:
I can load my kml file but I want to be able to define a marker for my points.  I can't find anywhere in the api which allows me to define the marker.

How can you define your own marker using a kml file with the esri javascript api?


texts: Setting CRS : TypeError: 'NoneType' object is not callable

texts:
I'm working with a naive geometry column imported from a csv file, like this:

      name       geometry                                           id
1     Abubilla  (POLYGON ((424405.897979997 4584795.70258297, ...   83
28    Abubilla  (POLYGON ((428405.900406701 4584295.67145889, ...   83
29    Abubilla  (POLYGON ((428405.919365087 4586795.67527349, ...   83
... ... ... ... ... ...
9808  Zarcero   (POLYGON ((424905.883573821 4582795.69571411, ...   40
9809  Zarcero   (POLYGON ((425405.880556747 4582295.69113956, ...   40
9810  Zarcero   (POLYGON ((425405.888148001 4583295.69267462, ...   40

When I try to set the CRS by the crs function:

gdf.crs({'init': 'epsg:25831'})

it yields this error message:

TypeError: 'NoneType' object is not callable

Could it be that Python doesn't recognize the geometry?
Is there another way to set the coordinates system for this "naive geometry"?