texts: Set default values when value is added to value table in Python Toolbox tool

texts:
Similar to Set default values for value table in Python Toolbox tool, I'm trying to populate a value table with default values. However for the linked question the value table is populated once an input feature class is input. I'm trying to add default values once the first value of the value table has been added. Here's my parameters:

class Tool(object):
    def __init__ (self):
        self.label = "Fast Field Calculate"
        self.description = ""
        self.canRunInBackground = False

    def getParameterInfo (self):
        inTab = arcpy.Parameter (
            displayName = "Input Table",
            name = "table",
            datatype = "DETable",
            parameterType = "Required",
            direction = "Input")

        calcs = arcpy.Parameter(
            displayName= "Update Field",
            name= "calcs",
            datatype= "GPValueTable",
            parameterType= "Required",
            direction= "Input")

        letters = string.ascii_uppercase [:numCol]
        calcs.parameterDependencies = ["table"]
        clmns = ["String", "Update Field"]
        clmns = [clmns] + [["String", "Calc Field '{}'".format (a)] for a in letters]
        clmns = clmns + [["String", "Logic Test"], ["String", "Equation"], ["String", "Null Values"]]
        calcs.columns = clmns
        calcs.filters [-1].list = ["TREAT AS NULL", "TREAT AS ZERO", "SKIP COMPUTATION"]

Here is the updateParameters function where I'm trying to create default values:

def updateParameters(self, parameters):
    inTab, calcs = parameters
    if inTab.altered and not inTab.hasBeenValidated:
        fc = inTab.valueAsText
        inFlds = arcpy.ListFields (fc)
        flds = ["-"] + [f.name for f in inFlds]
        for i in range (numCol + 1):
            if i == 0:
                fldFlds = list (flds)
                fldFlds.remove ("-")
                calcs.filters [i].list = fldFlds
                continue
    if inTab.altered and inTab.value and calcs.values:
        fc = inTab.valueAsText
        inFlds = arcpy.ListFields (fc)
        flds = ["-"] + [f.name for f in inFlds]
        for i in range (numCol + 1):
            if not calcs.values [-1] [i]:
                calcs.values [-1] [i] = "-" ##doesn't work - no default value displayed
            calcs.filters [i].list = flds
        if not calcs.values [-1] [-1]:
            calcs.values [-1] [-1] = "TREAT AS NULL" ##doesn't work - no default value

However, none of the values in the value table get populated. Is there a way to populate default values in a value table once the first value has been input?


texts: Some areas of label repeat, do not show, or appear offsets in Geoserver

texts:
I found that tile map of each region in geoserver2.9.2 of platform independent binary displays its name two or more times. As shown in the first picture.

I have done as How to avoid duplicate labels on each tile? [duplicate] and Polygon label repeated for each tile. However, some areas of label in sixth and seventh level do not show and some appear offsets. As shown in the second picture.

Here's my code:

<TextSymbolizer>
 <Geometry>
  <ogc:Function name="centroid">
     <ogc:PropertyName>the_geom</ogc:PropertyName>
  </ogc:Function>
 </Geometry>
</TextSymbolizer>

And then I have done as Labels on tiled map sometimes unavailable. It was same as above in seventh level.

Here's my code:

<VendorOption name="partials">true</VendorOption>

I have no idea how to modify it. Can anybody tell me?


texts: How to use Delaunay triangulation in QGIS?

texts: I have to create a Thiessen polygon for a site in New Jersey and all I have right now is the concentrations from samples we took in the area. What else do I need and how do I input it into QGIS to get my desired outcome? Which is to divide the total excavation into polygons based on concentrations? Any info will be a help.


texts: Dynamically add style to layer with user inputs

texts:
I am trying to create an SLD dynamically using parameters accepted from user. So how could I accomplish the task of adding the dynamically created SLD to my map layer?

Please tell me how to add style to layer dynamically, I am using GXP, OpenLayers, GeoEXT, GeoServer.


texts: QGIS: opening a shapefile without an exclusive file lock

texts:
UPDATE: I did some analysis of what happens with the locking and it looks like QGIS doesn't use an exclusive lock. It just requires (non-exclusive?) write lock to the file. So this C# code fails:

File.Open (fileName, FileMode.Open, FileAccess.Read, FileShare.Read);

the error being: System.IO.IOException : The process cannot access the file '...' because it is being used by another process

while this one goes through:

File.Open (fileName, FileMode.Open, FileAccess.Read, FileShare.ReadWrite);

So my question is basically off the mark and I'll move to close it.

In QGIS, when you open a shapefile in read-only mode, the shapefile still gets to be exclusively locked by QGIS, preventing other programs from reading it at the same time. This is pretty annoying and basically unnecessary, since you're not really editing the shapes or attributes of the file (and even if you were, the read access should still be allowed).

Is there a way to open the shapefile without locking it out for other software?


texts: Seeking geographic information for American and world parks and cities?

texts:
I am working on building out park data for my web site.  My site is about hiking, so I want to get a list of parks and cities in the world.

Where can I get/download/see such unified information?


texts: Lon/Lat coordinates for imported shapefile are displaying in the thousands

texts:
I'm having some trouble with my Italian regions shapefile.

Here I have imported my shapefile into QGIS but as you can see the Coordinate field at the bottom is showing very high numbers, and not the longitude and latitude that I am expecting. How can I see the correct longitude and latitude?

The file is downloaded from here http://www.istat.it/it/archivio/124086


texts: Selecting polygons within others that share a boundary

texts:
Thanks in advance for any help. I have a feature class in ArcGIS 10.2 with 455 polygons,they are municipal boundaries. We use it to determine permissions for users in another application and it's based on sharing a boundary. So, for example, Muni A touches Muni B and C, so Muni A's rights will be:

-Muni A

-All towns within

-Muni B

--All towns within

-Muni C

--All towns within

It's simple enough to select by location to determine the municipalities that share a boundary, however, I cannot seem to take that selection and do another based on those results to select the town polygons within B and C. Second to that, can you possibly script this to make it iterative so I can generate shp files for said permissions automatically?


texts: QGIS gives different resolution basemap in Print Layout

texts:
I use tiled background map, from online source with 512-pix tiles and try to use it as background for printed map, using Print Layout and PDF export features. Problem: it seems that the choice of tile zoom levels depends automatically from the DPI setting.

Basemap layer is defined simply as "XYZ tiles" layer with following URL: http://1.basemaps.cartocdn.com/rastertiles/voyager/{z}/{x}/{y}@2x.png

In the preview it uses screen dpi (70 perhaps), and I see this map on both map window and Print Layout. Tiles are here from zoom 16. And the street names are big and readable, as I need.

Now for for print I need 600 (or more) DPI PDF-s. When I set this as Export resolution of Layout and create PDF, then suddenly exporter takes tiles from zoom 18, so exported PDF has different map, as below. Probably QGIS tries to match print resolution with tile raw resolution, but this is bad idea for such raster tiles, as you end up with very small texts. This rescaling probably makes sense for satellites, but not with maps with scale-dependent styles like texts. The pattern is that with 300 DPI there is +1 zoom, with 600 +2 zoom etc.

In screenshots here the difference is not maybe so visible as it is in real printouts, but it is quite terrible.

Question: is there a way to adjust this multi-resolution raster layer resolution matching magic in QGIS Print Layouts?

This WMS question seems similar, maybe core issue is same: How to print a WMS map with chosen scale and WMS zoom level?

Btw, my quick remedy is now to use raster screenshot from Print Layout, but then I loose PDF vector with this, which would give much clearer vector texts and SVG vector icons on top of map.


texts: Polygons are unable to be edited and have no feature, although they are listed in attribute table with three features

texts:
I found that several lakes in the GRanD database (global reservoirs and dams) are unable to be selected for editing like normal polygons, and they have no feature at all, although they are shown in the attribute table with three features: GRanD ID, area and Poly_SRC.

May I ask:

Why are these polygons corrupted?
How to fix them to be normal polygons that can be edited?

Image 1: one lake polygon with no feature shown:

Image 2: the lakes shown in the attribute table with three features:

Image 3: The GRanD ID of each lake cannot be shown as label for each of these 25 mysterious polygons:

The shapefile is shared in this link.


texts: Reading JPEG2000 in R using gdal (installing gdal driver)

texts:
I currently write a code for automatic reading sentinel 2 data (JPEG2000 format), clipping it and save it a .tif. I use the language R and my system is Windows 10.

My problem is to install properly a gdal driver for JPEG2000. I tried to install the OpenJpeg library, but the installation is very poorly documented.

Can someone explain how to install the OpenJpeg library to use in R or has another suggestion how to read the JPEG2000 in R ?


texts: Incomplete geometry tool in QGIS 2.18 'Las Palmas'

texts: I just installed new version of QGIS which is 2.18 'Las Palmas'. Unfortunately, I cannot find vector conversion menu under geometry tool that I used to find in the previous version. I need to convert poly line to point, and polyline to polygon.


texts: Clipping tiff image with a bounding box with GDAL in Java

texts:
Since GDAL 2.1 GDAL and OGR utilities can be used as library functions, so this task is incredibly simple now:

from osgeo import gdal

ds = gdal.Open('original.tif')
ds = gdal.Translate('new.tif', ds, projWin = [-75.3, 5.5, -73.5, 3.7])
ds = None

Is there a similar library in Java to replace gdal_translate.exe call?


texts: Selecting features using an expression with PyQGIS

texts:
How can I select features with PyQGIS using an expression?

I tried to use an QgsExpression but the select method doesn't take it:

exp = QgsExpression("'ogc_fid' = 482")
cLayer = canvas.currentLayer()
cLayer.select(exp)

Is it possible and if so, how do I do it?


texts: Can variable be set to reference in_memory workspace object in ArcPy?

texts:
This is an offshoot question related to a question I posted earlier.  I have the following arcpy code that inputs a Python list of polygon geometry objects to the Union_analysis tool, which then writes the result to the in_memory workspace.  The result is then input to the Dissolve_management tool, which in turn writes the result to memory.  This all occurs near the end of a function, and the final goal here is to return the single geometry from the Dissolve to the calling function, where it can be passed along.  Hence my question:  How can (can I?) I set the Dissolve output to a variable so that I can return it as a geometry object?  The first thing that will be done with it is to read its geometry into a point array to use as the inner ring for another polygon.

    if len(lstHSGeoms) > 1:
        #Union geometries into one geometry
        arcpy.Union_analysis(lstHSGeoms, "in_memory\\lstHSUnionGeom", "ONLY_FID")
        arcpy.Dissolve_management("in_memory\\lstHSUnionGeom","in_memory\\lstHSOneGeom", None, None, "SINGLE_PART")
        hsOneGeom = lstHSOneGeom

The last line is the problem here.  No surprise, arcpy complains that the global name lstHSOneGeom is not defined.  Putting in the whole in_memory path as the assignment value, i.e. "in_memory\\lstHSOneGeom" just assigns the simple string to the hsOneGeom variable.  There is no reference to the in_memory object.

I originally coded this block as the ArcGIS help files would have you do, which is to create empty arcpy Geometry objects and to set them as the outputs of the tools.  When then used as inputs, they were rejected as invalid.  At least I could keep a reference to the geometry objects in a variable, but the tools just wouldn't recognize them as valid input.  So, I went with in_memory workspace update that you see in the snippet above.  (Thanks to @dmahr for this solution in the previous question referenced above.)  I now don't get errors from the tool, but I do get an error from the last line.  Any suggestions?

I'm using ArcGIS Desktop (ArcView) 10.0, Windows 7, 8 GB RAM.


texts: How to extract NetCDF raster values to time enabled shapefiles?

texts:
Does anyone know a good way to extract NetCDF raster values to each data point from a time enabled shapefile at the correct time?

I have tried using a python script tool and a loop that steps through time but it seems that the loop only works for the shapefile and not for the NetCDF raster - all shapefile values (from all time step intervals) end up getting assigned raster layer values just from the first time step interval. Both the raster and the shapefile are time enabled with the same time step interval and I have checked with Time Slider that the values of both change through time correctly.

Here is what I have tried:

mxd = arcpy.mapping.MapDocument('SST.mxd')
df = arcpy.mapping.ListDataFrames(mxd)[0]
rasterLayerName = "SST"
rasterLayer = arcpy.mapping.ListLayers(mxd, rasterLayerName, df)[0]
plasticLayerName = "Plastic"
plasticLayer = arcpy.mapping.ListLayers(mxd, plasticLayerName, df)[0]
df.time.currentTime = df.time.startTime
while df.time.currentTime <= df.time.endTime:
   uniqueName = "Plastic_" + df.time.currentTime.strftime("%Y_%m")
   outPointFeatures = "U:/gis/" + uniqueName
   arcpy.sa.ExtractValuesToPoints(plasticLayer, rasterLayer,
         outPointFeatures, "NONE", "VALUE_ONLY")
   df.time.currentTime += df.time.timeStepInterval


texts: Adding Geographic Coordinate System and Projected Coordinate System layers

texts: So I have two layers that I am adding into ArcMap.  I know for a fact that when you added the first layer with Projected Coordinated System, the 2nd layer, which has a geographic coordinate system will project on the fly. My question is what would happen if the Geogrpahic coordinate System is added first. What would happen and what are the necessary steps to have them all fit properly in the same coordinate system ?


texts: ArcPy: Initialize zero valued grid Raster and iteratively increment

texts:
I'm jumping right into ArcPy to help out a friend who needs to automate some ArcGIS work. The problem at hand can be reduced to the following pseudocode:

accumulationRaster = SomeZeroedRaster
for each Raster pair, R, S, pertaining to rain and snow:
   accumulationRaster += Con(S == 0, Con(R >= 1, 1, 0), 0)

In other words, iterate through all data of rain and snow of each day. Perform a logical operation on each cell pair, so that the corresponding cell in the accumulationRaster gets incremented if no snow (S == 0) and some rain (R >= 1).

Having no luck so far with various attempts, I thought I'd ask here, as the problem should have an easy solution. Full working code with the missing gaps pointed out below. Most of it is just trivial file path set up, and the above pseudocode should suffice:

The missing pieces are numbered 1, 2, and 3.

How to initializase a zero raster w/proper dimensions.
How to make a copy of this zero raster. If 1 is simple, just do this instead.
Increment accumulation raster.

import os
import arcpy
from arcpy import env
from arcpy.sa import *

in_workspace = "C:/Data/ProjectFolder"
os.chdir(in_workspace)
arcpy.env.workspace = in_workspace

maaneder = ["01_januar", "02_februar", "03_mars", "04_april",
            "05_mai", "06_juni", "07_juli", "08_august",
            "09_september", "10_oktober", "11_november", "12_desember"]

################ 1 ###################
baseAccum = #1: ?? (zero Raster w/correct type & dimensions)
######################################

outputPath = "rainDaysInGrowthSeason/"
for year in range(1997,2015):
    for month in range(1, 13):
        ################ 2 ###################
        monthAccum = #2: ?? (Copy of baseSum)
        ######################################

        for day in range(1,32):
            # e.g 1997/09_september/
            fpath = str(year) + "/" + maaneder[month-1] + "/"
            # e.g. 1997_09_08.asc
            filesuffix = str(year) + "_" + str(month).zfill(2) + "_" + str(day).zfill(2) + ".asc"
            # e.g. Snoedybde/1997/09_september/sd_1997_09_08.asc
            snoefile = "Snoedybde/" + fpath + "sd_" + filesuffix
            # e.g. rrl/1997/09_september/rrl_1997_09_08.asc
            rrlfile = "rrl/" + fpath + "rrl_" + filesuffix

            hasSnoe = os.path.isfile(snoefile)
            hasRrl = os.path.isfile(rrlfile)
            if (not hasSnoe or not hasRrl):
                continue

            snoeRaster = Raster(snoefile)
            rrlRaster = Raster(rrlfile)

            ############## 3 #####################
            monthAccum = #3: ??
                         # something a'la:
                         # += Con(snoeRaster == 0, Con(rrlRaster >= 1, 1, 0), 0)
            ######################################

        # e.g. rainDaysInGrowthSeason/1997
        outputDir = outputPath + str(year)
        if (not os.path.isdir(outputDir)):
            os.mkdir(outputDir)
        monthAccum.save(outputDir + "/r_" + str(year) + "_" + str(month).zfill(2))


texts:
I  would like to convert an img to grid in ArcMap 10.2. I tried it in two ways:

In Conversion tools with the Raster to other format tool, where I adjusted the GRID format and the output workspace was a gdb
In ArcCatalog with right click on the image I choosed Export --> Raster to Different Format. Here I didn't adjust any extension and I tried to save the output raster into a gdb and also into a simple folder.

My problem is that the result is in FGDBR format in every case. What could be the problem?


texts: Anyone know of a good reference (either pay or open-web) that details recommended country projections?


texts: Unexpected records added after spatial join

texts:
In the following spatial join, I want to join pd_subdistricts to parcel and keep both matching and unmatching records.  Since parcel has 292,109 records, the output should have the exact same amount, given that geometries are valid.  After making the geometries valid, I get an unexpected result.   The output - test has 249 more records than parcel  How is this possible?

I have performed the exact same query to join other tables and the resulting number of records is always equal to parcel

Input   parcel          has 292,109 records
Input   pd_subdistricts has 1,083 records
Output  test            has 292,358 records

Using this query,

CREATE TABLE test AS
SELECT t.*
FROM parcel AS t
LEFT JOIN pd_subdistricts AS m
  ON ST_Within(t.geom, m.geom)


texts: Geoservers communicating with other Geoservers

texts: I am able to connect to multiple Geoservers using the Geonetwork but is it possible for two Geoservers to communicate with each other without the Geonetwork?


texts: PostGIS 2.0 ST_Transform with Rasters cannot find custom SRID

texts:
I'm trying to calculate an average value in a buffer around a point.

Given a raster in EPSG 4326 and a similar point, this is fine, except that the values will be skewed as they are not projected in an equal-areas projection.

So - I've added to my spatial_ref_sys the SRID '97099' from saptialreference.org:
http://spatialreference.org/ref/sr-org/7099/

This works fine for transforming points [ edit - added example ]:

-- get all the diff areas of diff pix values in one tile in which a point lies
SELECT ST_Area((ST_DumpAsPolygons(r.rast)).geom),
       (ST_DumpAsPolygons(r.rast)).val,
       r.rid
FROM rasters AS r
INNER JOIN (
    SELECT ST_Transform(ST_SetSRID(ST_MakePoint(-2, 51), 4326), 97099) as pt
) AS p
ON ST_Intersects(r.rast, p.pt) AND r.filename = 'someraster.tif';

But using it with the raster fails with an error saying it cannot find SRID 97099 in spatial_ref_sys: [ edit - added error message ]

ERROR: rtpg_getSR: Cannot find SRID (97099) in spatial_ref_sys
SQL state: XX000
Context: SQL function "st_transform" statement 1

A SELECT ... WHERE SRID = 97099 confirms it is there... So I'm a bit lost. Is it querying from somewhere else, or maybe ignoring anything that isn't EPSG? I've followed the example here, for rasters:
http://www.postgis.org/documentation/manual-svn/RT_ST_Transform.html

Here's the full test SQL, the commented out lines show what I want to do, and the line above them will actually work (but is wrong, as it will still be 4326)

DROP TABLE IF EXISTS tmp_points;
CREATE TEMP TABLE tmp_points ("gid" int, "point" geometry);
INSERT INTO tmp_points ("gid", "point") VALUES
  (100231, ST_MakePoint(172.3, -43.76667)),
  (100233, ST_MakePoint(171.71667, -43.91667)),
  (100246, ST_MakePoint(170.92009, -44.09834)),
  (100290, ST_MakePoint(170.95009, -44.44835)),
  (100497, ST_MakePoint(173.28271, -42.77409)),
  (100498, ST_MakePoint(171.90983, -43.9236)),
  (1, ST_MakePoint(-2.136154, 49.220288))
;

-- get the average value arount a point for a given radius
--
SELECT
  gid,
  filename,
  lon,
  lat,
  CAST(SUM(ST_Area((foo.gv).geom)*(foo.gv).val)/SUM(ST_Area((foo.gv).geom)) AS decimal(9,7)) AS avgimr,
  COUNT(foo.gv) as parts
FROM (
  SELECT
  ST_X(p.point)::NUMERIC(9, 5) AS lon,
  ST_Y(p.point)::NUMERIC(9, 5) AS lat,
  p.gid,
  sn.filename,
  ST_Intersection(
    sn.rast,
    -- THIS COMMENTED LINE BELOW FAILS, THE ABOVE WORKS (is wrong)
    -- ST_Transform(sn.rast, 97099, 'Bilinear'),
    ST_Buffer(ST_Transform(ST_SetSRID(p.point, 4326), 97099), 25000)
  ) AS gv
  FROM rasters sn, tmp_points p
  WHERE ST_Intersects(
    ST_Buffer(ST_Transform(ST_SetSRID(p.point, 4326), 97099), 25000),
    sn.rast
    -- THIS COMMENTED LINE BELOW FAILS, THE ABOVE WORKS (is wrong)
    -- ST_Transform(sn.rast, 97099, 'Bilinear')
  )
  AND sn.filename = 'someraster.tif'
) AS foo
WHERE (foo.gv).val >= 0 AND (foo.gv).val <= 10
GROUP BY filename, gid, lon, lat
ORDER BY filename, gid, lon, lat;

Many thanks


texts: Can a geojson multigeometry contain a multipoint among other geometries?

texts:
Can a Geojson MultiGeometry contain a Multipoint array of points, or is it necessary to convert them into individual Points before adding to the geometry collection?

This related question suggests it is allowed.
Can a GeoJSON GeometryCollection contain another Collection?


texts: Convert WKT coordinate reference system to PROJ4

texts:
I'm currently working on a better support of LAS 1.4 file format in the rlas R package. In LAS 1.4 format the CRS of the point cloud can be stored as a WKT string. For example this is what I read from a LAS 1.4 file (raw data).

"COMPD_CS[\"Projected\", PROJCS[\"UTM_10N\", GEOGCS [ \"WGS84\", DATUM [ \"WGS84\", SPHEROID [\"WGS 84\", 6378137.000, 298.257223563 ], TOWGS84 [ 0.000, 0.000, 0.000, 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000 ] ], PRIMEM [ \"Greenwich\", 0.000000 ], UNIT [ \"metres\", 1.00000000] ], PROJECTION[\"Transverse_Mercator\"], PARAMETER[\"Latitude_of_Origin\",0.0000000000], PARAMETER[\"Central_Meridian\",-123.0000000000], PARAMETER[\"Scale_Factor\",0.9996000000], PARAMETER[\"False_Easting\",500000.000], PARAMETER[\"False_Northing\",0.000], UNIT [ \"metres\", 1.00000000]] ], VERT_CS[\"NAVD88 (Geoid03) ContUS\", VERT_DATUM[\"./Resources/CoordSysData/navd88_geo03_contus.bin\", 1 ], UNIT [ \"metres\", 1.00000000] ] ]"

How can I interpret this string in R? In the R spatial ecosystem all packages rely on PROJ4. How I can convert this WKT string into PROJ4 string?

geos::readWKT does not understand this string.


texts: Shapefile of South Tyneside borough (north-east England) boundary

texts: I'm looking for a free-to-download shapefile (polygons) of South Tyneside borough. South Tyneside is a metropolitan borough in NE England. Can anyone help?


texts: Is there a limit to the number of GCPs that can be used?

texts:
I am trying to georeference an image using many (>1,000,000) GCPs (derived by matching two images using triangulation). However, after many attempts doing this in different ways I think I found that there is a limit to the number of GCPs that I can use.

ds = gdal.Open("tmp_GCP.TIF", gdal.GA_Update)
wkt = ds.GetProjection()
ds.SetGCPs(gcp_list_gdal[:10923], wkt)
print(ds.GetGCPCount()) # output 10923
ds.FlushCache()
ds = None
gdal.Open("tmp_GCP.TIF").GetGCPCount() # output 0

When I set 10923 GCPs I obtain a file that seemingly has zero GCPs and when I set this to 10922 I do get a file with 10922 GCPs. This seems to be independent of the width/height of the raster.

Also, I tried to set the GDAL cache to a higher value, but this also didn't have an effect.

gdal.SetConfigOption('GDAL_CACHEMAX', '1000000')

Edit:
Just found out that I am getting the following error in the terminal when running this code in a Jupyter Notebook (which I don't get when I use gcp_list_gdal[:10922]).

Warning 1: TIFFFetchNormalTag:Incorrect count for "GeoTiePoints"; tag ignored
Warning 1: TIFFFetchNormalTag:Incorrect count for "GeoTiePoints"; tag ignored

I have tried this on Ubuntu 16 and Ubuntu 18, using Gdal 2.2, 2.4 and 3.0, but I always get the same error.

Edit 2: Because of the error I thought it might be a problem with the amount of GCP information that can be stored in a GeoTIFF file, so I tried another approach in python (trying not to store the DS before applying a warping), but it seems the same problem occurs (since it does produce an output using 10922 GCPs but not when using 10923 GCPs).

ds = gdal.Open('tmp_GCP.TIF', gdal.GA_Update)
gcp_srs = osr.SpatialReference()
gcp_srs.ImportFromEPSG(4326)
wkt = gcp_srs.ExportToWkt()
ds.SetGCPs(gcp_list_gdal[:10922], wkt)
dstDS = gdal.Warp('TEST.TIF', ds, format='GTiff', tps=True, dstNodata=0,
                  srcNodata=0, resampleAlg=gdal.GRA_Bilinear)

Edit 3: I found out about the option to use geolocation arrays https://gdal.org/development/rfc/rfc4_geolocate.html which seems very interesting. As I understand it you can basically create two bands of the same size as the image and indicate for each pixel the desired longitude and latitude of the pixel. This should work better with a large number of GCPs so I am going to try this.


texts: Bringing in UTM NAD17 data from excel into CartoDB

texts:
Any ideas how to make a map from an excel file that has UTM Zone 17 NAD27 x,y's?  This is what I wrote in my SQL:

SELECT * FROM storklocations2
UPDATE storklocations2 SET the_geom = ST_Transform(ST_SetSRID(ST_Point(easting, northing),2031),4326)

It gives me an error.


texts: What is a correct overpass-turbo query for getting all streets in a city?

texts:
I want to get all streets in NYC using http://overpass-turbo.eu/. I tried this:

[out:json]; area[name = "New York"]; (node(area)[highway=street]; ); out;

However it returns

{
  "version": 0.6,
  "generator": "Overpass API 0.7.55.1009 5e627b63",
  "osm3s": {
    "timestamp_osm_base": "2019-11-13T19:26:03Z",
    "timestamp_areas_base": "2019-11-13T18:05:02Z",
    "copyright": "The data included in this document is from www.openstreetmap.org. The data is made available under ODbL."
  },
  "elements": [

  ]
}

No elements. However this query:

[out:json]; area[name = "New York"]; ( node(area)[amenity=cinema]; node(area)[highway=street]; ); out;

for getting streets and cinemas, works:

{
  "version": 0.6,
  "generator": "Overpass API 0.7.55.1009 5e627b63",
  "osm3s": {
    "timestamp_osm_base": "2019-11-13T19:29:02Z",
    "timestamp_areas_base": "2019-11-13T18:05:02Z",
    "copyright": "The data included in this document is from www.openstreetmap.org. The data is made available under ODbL."
  },
  "elements": [

{
  "type": "node",
  "id": 344994897,
  "lat": 41.7680892,
  "lon": -73.9291000,
  "tags": {
    "amenity": "cinema",
    "created_by": "Potlatch 0.10f",
    "name": "Roosevelt Theater"
  }
},
{
  "type": "node",
  "id": 344996199,
  "lat": 41.7727941,
  "lon": -73.9292159,
  "tags": {
    "amenity": "cinema",
    "created_by": "Potlatch 0.10f",
    "name": "Hyde Park Drive-In"
  }
},
{
  "type": "node",
  "id": 368043485,
  "lat": 40.7259416,
  "lon": -73.9973321,
  "tags": {
    "addr:state": "NY",
    "amenity": "cinema",
    "ele": "10",
    "gnis:county_name": "New York",
    "gnis:feature_id": "2082840",
    "gnis:import_uuid": "57871b70-0100-4405-bb30-88b2e001a944",
    "gnis:reviewed": "no",
    "name": "Angelika Film Center",
    "source": "USGS Geonames",
    "wikidata": "Q45104731"
  }
},
...

How should I modify the initial query to get the streets?

EDIT:

I tried

[out:json];
area[name = "New Hampshire"]->.a;          // Redirect result to ".a"
out body qt;
(
  way
    (area.a)                    // Use result from ".a"
    ["highway"="living_street"];
);
out body qt;
>;
out skel qt;

but it didn't select all of the streets, just several. And how should I specify that I mean a state by New Hampshire? I need all streets with postcodes (zip).

this doesn't display any streets at all:

[out:json];
area[name = "Acworth"]->.a;          // Redirect result to ".a"
out body qt;
(
  way
    (area.a)                    // Use result from ".a"
    ["highway"="living_street"];
);
out body qt;
>;
out skel qt;

Although it should at least display these: https://geographic.org/streetview/usa/nh/sullivan/acworth.html


texts: Merging Urban Paths and ITN data in a network

texts: I have an ITN network created within ArcMap. I also have urban paths data for the same area. I would like to merge both datasets into one network. I would like the network to be multimodal. I cannot figure out how to do this. I have tried using productivity suite, but I cannot create a network using both Urban Paths and ITN simultaneously.


texts:
I am trying to interpolate a point layer holding data of type double, I have gone back to older projects where I have successfully interpolated before and these don't work either (resulting in nan).

I have checked and the layers are of the same CRS, what doesn't make sense is the older project also not interpolating.

Is anyone else recently experiencing the same issue, if so have you managed to get around it?


texts:
I am attempting to make a simple select and zoom python script for my address points. Currently the script will complete successfully but will not select or zoom to the entered point. I am using two parameters that have been set in the Script Tool. Address_Number (Long) and Address_Name (String) which both have identical spellings/data types in Properties>>Parameters

Here's the code:

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
addLayer = arcpy.mapping.Layer("G:\Layers\GIS_SDE.AddressPoints.lyr")
arcpy.mapping.AddLayer(df, addLayer, "TOP")
Address_Number = arcpy.GetParameterAsText(0)
Address_Name = arcpy.GetParameterAsText(1
arcpy.SelectLayerByAttribute_management("GIS_SDE.AddressPoints", "NEW_SELECTION", "A_NBR = %Address_Number% AND ADDR = %Address_Name%")
df.zoomToSelectedFeatures()


texts:
I am using QGIS 3.14. I want to organize my columns in the attribute table but I cant change the virtual fields columns place.
It is always last. Can I put it between two non-virtual fields columns?


texts: Distinguishing two streets with the same name?

texts:
I am trying to get road data using the Overpass API. I want a user to be able to find a street using its name then retrieve the data for the nodes to draw a line on the map. This is all fine using:

[out:json][timeout:25];
(
  way["name"="<streetName>",i](51.388370, -2.421840, 51.410317, -2.361116);
);
(._;>;);
out body;

The problem is when a road has the same name as another in the bounding box. For example in the UK "High Street" is a common street name. I can't see any tags (even using verbose output) that give me more information about where the street is. Obviously there is the actual geographic node data but this isn't human-friendly.

Is there any way to work out any more information about a way (specifically a street/highway) to distinguish them to humans?


texts:
I am new to leaflet & its geojson features.

I am trying to test leaflet's functionalities and have built a map and its layer:

var map = new L.map('map', {scrollWheelZoom:true, doubleClickZoom:true, boxZoom:true, zoomControl:true});
var osmTiles = new L.TileLayer('http://localhost/googleapi/images/image1/tiles/{z}_{x}_{y}.png', {
maxZoom: 7,     minZoom: 3
}).addTo(map);

Now i want to draw a second layer (geoJson) upon it with a marker, so following the tutorila (cloudmate):

var geojsonLayer = new L.GeoJSON();
var geojsonFeature = {
"type": "Feature",
"properties": {
    "name": "Coors Field",
    "amenity": "Baseball Stadium",
    "popupContent": "This is where the Rockies play!"
},
"geometry": {
    "type": "Point",
    "coordinates": [83.06877, -156.35742]
}
};
geojsonLayer.addData(geojsonFeature);

I am testing in my local PC (have not installed node + built leaflet) so maybe i am missing something.


texts: Setting arcpy.env.workspace to ArcSDE connection in separate module on ArcGIS for Server?

texts:
I have problem defining arcpy.env.workspace in a separate module when publishing the script as a Geoprocessing Service on arcgis server

For instance:

my main tool script:

from myModule import myTool

workspace = "Database Connections\\GEOP011@demma18.sde"
tableName = "DEMMA18.MyTable"

myTool.execute(workspace, tableName)

myModule.myTool:

import arcpy

def execute(workspace, tableName)
    arcpy.env.workspace = workspace
    doSomething(tableName)

def doSomething(tableName)
     with arcpy.da.SearchCursor(tableName, ["OBJECTID"]) as cursor:
            for row in cursor:
                  print(row)

When published to the server, this script will return error : cannot open DEMMA18.MyTable

However, if I put the same code in a single .py file it runs on the server without error. Any Idea what is going wrong?


texts:
I'm trying to make a hex grid net with the shape of polygon. Unfortunately, the grid net is in the rectangle shape (see the screenshot).

Any advice how to do it?


texts: Using Geocoding and creating the right locator in ArcGIS for Desktop?

texts:
I have a set of addresses that I am wanting to geocode. I have available a parcel/points layer consisting of addresses in the region which includes the fields - street address number, road name and road type along with a full string of these in one field.

e.g. Address number - 12    Road name - Hunterway Road   Road type - Road     Full address - 12 Hunterway Road, Chartwell,                          Hamilton

I am hoping to use this layer above to create an address locator with the locator style - US Address - Single House.

USing this created locator I would like to geocode the address I have available which are currently in a CSV file format and have the following fields:

Add1  -        12 Hunterway Road
Add2     - Chartwell
City_add - Hamilton
Address_string - 12 Hunterway Road, Chartwell, Hamilton, 3210

Currently I am not getting this to work and was wondering if the created locator is the wrong type or whether I have to do some table edits with either the data used to create the address locator or with the data I am hoping to geocode.


texts:
Is it possible to set the delay timer when hovering the cursor over maptips?

It seems like there's a one second delay as shown below but I would like to have it appear instantly instead.

Using QGIS 2.18.17 and 3.0.1.


texts: Extracting pixels in multiple rasters based on the values of another raster's pixels for time series analysis

texts:
I'm working with 52 rasters of a specific region over a period of one year in python. I also have a mask raster of the same area that is a binary color (250 black pixels for disparate areas of interest, nothing for everything else).

I'm trying to use a black pixel from the mask raster as bounds to grab a single pixel over time (52 rasters/weeks) as a numpy array, and I'd like to do this for all 250 black pixels, grabbing their corresponding pixels over time into separate arrays.

Unfortunately I cannot seem to find a method to pull this off. The closest thing I've found is using a shapefile to clip a raster. Any suggestions on how to achieve this in python would be golden.


texts: Snippets applicable for ArcObjects 10.2.1

texts:
I am trying to create ArcGIS desktop add-in using Visual Studio Professional 2012. But when I use ArcGIS Snippets Finder to insert code, I found the applicable versions as follows:

' Applicable ArcGIS Product Versions:
' 9.2
' 9.3
' 9.3.1
' 10.0

But I am using ArcGIS 10.2.1, and the ArcObjects obtained from ArcGIS 10.2.1 Desktop installation CD.

I want to know where can get the snippets applicable for ArcGIS 10.2.1 (or above). And what version of Visual Studio can support ArcObjects for 10.2.1?


texts: ST_Intersects performance issue and not using spatial index with simple update query

texts:
I am having a performance issue with ST_Intersects. I know similar questions have been asked several times already, but I cannot seem to make use of the past answers in my case. I am trying to perform a very straightforward update to an addresses table using St_Intersects. I have 500,000 addresses (points) and I want to capture the ID of the flood zones which each address intersects. There are 5,000 flood zones (polygons).
This is my query:
UPDATE nfm.nfm_addresses a
SET fld_zone_id = f.ogc_fid
FROM nfm.s_fld_haz_ar_36103 f
WHERE ST_Intersects(a.geom, f.geom);

I cannot index the flood zone geometry due to the size of some polygons exceeding index-able limits. I try anyway and let the query run overnight to no avail. Here is the EXPLAIN:
"Update on nfm_addresses a  (cost=0.00..872514959.88 rows=1097483 width=184)"
"  ->  Nested Loop  (cost=0.00..872514959.88 rows=1097483 width=184)"
"        Join Filter: ((a.geom && f.geom) AND _st_intersects(a.geom, f.geom))"
"        ->  Seq Scan on nfm_addresses a  (cost=0.00..14264.08 rows=587308 width=174)"
"        ->  Materialize  (cost=0.00..1818.09 rows=5606 width=42)"
"              ->  Seq Scan on s_fld_haz_ar_36103 f  (cost=0.00..1790.06 rows=5606 width=42)"

I then used St_Subdivide to get the flood zones small enough to be indexed. This resulted in 50,000 flood zone polygons instead of the original 5,000. I am able to index the geometry in this case and attempt the same query. Here is the EXPLAIN:
"Update on nfm_addresses a  (cost=0.00..7579370470.20 rows=9533770 width=184)"
"  ->  Nested Loop  (cost=0.00..7579370470.20 rows=9533770 width=184)"
"        Join Filter: ((a.geom && f.geom) AND _st_intersects(a.geom, f.geom))"
"        ->  Seq Scan on nfm_addresses a  (cost=0.00..14264.08 rows=587308 width=174)"
"        ->  Materialize  (cost=0.00..8570.49 rows=48699 width=42)"
"              ->  Seq Scan on s_fld_haz_ar_36103 f  (cost=0.00..8326.99 rows=48699 width=42)"

No geometry index appears to be used. I've tried rearranging my query in several ways using CTEs and/or subqueries to no avail. What am I missing?
I'm using PostgreSQL 11.3, compiled by Visual C++ build 1914, 64-bit. I'm using POSTGIS="2.5.2 r17328" [EXTENSION] PGSQL="110" GEOS="3.8.0-CAPI-1.13.1 " PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.4.4, released 2020/01/08" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER.
My machine is a top of the line desktop with an i9, 64GB RAM, etc.


texts:
I'm using MapInfo 11.0.4 Release Build 404

Whenever I want to save a workspace file in the MapInfo workspace format (*.mws) it saves the file as expected and with the right file extension. However, when I look at the content of the file in a text viewer, I can see that it is a .wor file in reality. It's not saving with the XML format that should be in a *.mws file.

Since there is no other way to save workspace I cannot create the MWS format at all and I need it to use with the MapInfo Spacial Server product.

Any help would be appreciated.


texts: Hosted tile layer's pop-up won't work in AGOL web map

texts: Upon opening one of my web maps I realized that the tile layers used in this map no longer display their pop-ups. These pop-ups were given to the tile layers from feature (polygon) layers. When trying to configure the tile layers' pop-ups all appears fine. This was not an issue till recently, and the web map operated as intended. Has a recent update caused this bug, or could this be something on my side?


texts: How to get GeoJSON property name in Google Maps JavaScript API?

texts:
I would like to show the tooltip with the unknown property name and its value in Google Maps JavaScript API application. GeoJSON will be uploaded by the user and that's why property names are unknown. I have JS code like this:
map.data.addListener('click', function (event) {
      var value_name1 = ?
      var value_1 = event.feature.getProperty(PutHereDynamicValue);

      let html = value_name1 ":" value_1;
      infowindow.setContent();
      infowindow.setPosition(event.latLng);
      infowindow.setOptions({ pixelOffset: new google.maps.Size(0, -10) });
      infowindow.open(map);
    });
  }

How can you get dynamically from GeoJSON its property name and the value? Should I parse GeoJSON and extract properties keys?


texts: Intensity Wise Gradient is not working in leaflet.heat

texts:
I am using leaflet.heat for creating heatmap with leaflet API.I have a highest intensity of 1.0 and lowest 0.0. When I am drawing a heatmap in the highest zoom its all are showing the same color. I have plotted the value over the heat color. Here you can see what I am getting in the attached Image. Here is the code I am using. Seems like Intensity is not working and using all value as same.

var heatmap = L.heatLayer(allHeats, {
                radius: 20,
                max: 1.0,
                blur: 15,
                gradient: {
                    0.0: 'green',
                    0.5: 'yellow',
                    1.0: 'red'
                },
                minOpacity: 0.7
            }).addTo(map);

allHeats Array data

[23.703161, 90.414101, 1],
[23.722657, 90.412525, 0.6],
[23.701936, 90.415782, 0.5],
[23.722099, 90.413981, 0.9],
[23.703582, 90.414536, 0.3],
[23.699689, 90.421493, 0.1],
[23.705727, 90.412281, 0.8],
[23.720692, 90.414056, 0.2],
[23.706133, 90.412448, 0.8],
[23.704677, 90.41166, 0.7],
[23.703912, 90.413993, 0.4]


texts: How to display multiple attributes in a QGIS Composer legend?

texts:
I have a polygon layer which has a number attribute and a text attribute.
When I make a legend in the composer, it automatically chooses the number to display while I would like to have both the number and the text.

How do I display two attributes in a legend? If not possible, how do I select which attribute is displayed in the legend?

I am not looking to have a multi-attribute label in the main window, I wish to have a multi-attribute label in the legend of the composer, or choose which attribute is displayed in the legend.

QGIS 2.8 Win8 64 bit.


texts: Display labels for all point features in ArcGIS

texts: I might be missing something silly, but is there a way to guarantee that each point label in my point layer gets displayed?


texts: Move to next iteration in event of empty feature class

texts:
I'm working on a model that processes lidar data from a public layer file, one quadrangle at a time, as I don't have a powerful enough machine, or enough storage to do it all at once. I am looking for cliff bands based on slope, through raster calculator I convert the slope layer into a binary representation, showing only areas above a certain steepness threshold, next I contour the binary slope representation, buffer the contours, then am performing zonal statistics on the buffered contours. The issue I'm running into, is that if there are no cliffs I get an error at the zonal statistics, because of an empty feature class.

Is there a way to make my model move to the next iteration when a particular error arises, so that I may process my whole area of interest, despite some of the individual chunks I analyze have no cliffs/will produce an error?

Dissolving...
WARNING 000953: Extent of the overlay operation will be empty.
WARNING 000117: Warning empty output generated.
Succeeded at Fri Jun 22 14:27:16 2018 (Elapsed Time: 0.66 seconds)
Executing (Zonal Statistics): ZonalStatistics "C:\Users\NAME\Desktop
\CLIFFFINDERPURGEFOLDER\CLIFFFINDERPURGEGDB1.gdb\TIGHTANALYSISBUFFER_2"
OBJECTID "C:\Users\NAME\Desktop\CLIFFFINDERPURGEFOLDER
\CLIFFFINDERPURGEGDB1.gdb\DEMCLP2" "C:\Users\NAME\Desktop
\CLIFFFINDERKEEPFOLDER\CLIFFFINDERKEEPERGDB1.gdb\MAXCLIFFHEIGHT_2" RANGE
DATA
Start Time: Fri Jun 22 14:27:16 2018
ERROR 999999: Error executing function.
The table was not found. [VAT_DEMCLP2]
ERROR 010151: No features found in C:\Users\NAME\Desktop
\CLIFFFINDERPURGEFOLDER\CLIFFFINDERPURGEGDB1.gdb\TIGHTANALYSISBUFFER_2.
Possible empty feature class.
Failed to execute (Zonal Statistics).


texts: Strange behavior when writing coordinates to an ESRI shapefile

texts:
I am using python code straight out of the GDAL/OGR Cookbook to write lat/lon coordinates (stored in a numpy array of tuples) to an ESRI shapefile like this:

#using almost entirely GDAL cookbook code here,
#but have already done my affine transform

multipoint = ogr.Geometry(ogr.wkbMultiPoint)
for i in range(len(coordtrees)):
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(coordtrees[i,0], coordtrees[i,1])
    multipoint.AddGeometry(point)

shpDriver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists('points.shp'):
    shpDriver.DeleteDataSource('points.shp')
outDataSource = shpDriver.CreateDataSource('points.shp')
outLayer = outDataSource.CreateLayer('points.shp', geom_type=ogr.wkbMultiPoint)
featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(ogr.Geometry(ogr.wkbMultiPoint))
outLayer.CreateFeature(outFeature)
outFeature = None

When I attempt to check the features (also pulled from the GDAL/OGR cookbook):

file = "/pathname/points.shp"

driver = ogr.GetDriverByName('ESRI Shapefile')

dataSource = driver.Open(file, 0)

if dataSource is None:
    print(f'Could not open {file}')
else:
    print(f'Opened {file}')
    layer = dataSource.GetLayer()
    featureCount = layer.GetFeatureCount()
    print(f"Number of features in {os.path.basename(file)}: {featureCount}")

I get back that I have 0 features in the .shp file (same when I view it in QGIS). However, when I rerun the first script, ceteris paribus, and check my features, I have my Multipoint file with one feature (that also displays in QGIS)! Why is this happening?


texts: where is the v.transects algorithm in qgis?

texts:
I am interested in using the v.transects GRASS command. I am however not yet proficient with using GRASS or python in QGIS environment so I would rather do it the old-fashioned way (i.e using the GUI).

When I open my geoprocessing toolbox in QGIS however I don't find it among the listed GRASS functionality. Why is this particular command missing? How can I add it? If not will I have to resort to doing this operation through a command line of some sort (prefer to avoid)

I did end up trying to open the layer in GRASS and running the following command:

v.transects --help

To which I got the following error:

  'v.transects' is not recognized as an internal or external command,
  operable program or batch file.

I know the session environment is OK because when I try to run another command, such as v.buffer --help, I get valid output...

Does the v.transects command even exist? if not why does it appear in the documentation?


texts: How to convert OSM Data to PostgreSQL for Routing Purposes?

texts:
I'm creating a routing application using OpenStreetMap. Currently, my task is to convert the data from an OSM file to PostgreSQL.

I read about osm2pgsql, but the problem is it says there that osm2pgsql is not suited for routing. I read about Osmosis too, but there is something about rails port, which I think means I need to learn Ruby in Rails. I tried studying RoR in the past and I found it too difficult for me to learn.

I'm using PHP for my application. How can I able to extract the data from the OSM file and export it to my PostgreSQL. I know I might be missing something here, but I'm a newbie when it comes to mapping application. I hope someone can help me. Thanks! :D


texts: How to edit the attribute tabel of a shapefile when the enable editing is greyed out?

texts:
I am not really confident with QGIS.

I am looking for at way to  edit the attribute table of a shapefile, but the enable editing button is greyed out.

How do i fix it?

My end goal is to add a column to the file I have (which is a world map), so I can color code the world's nations according to their CO2e emissions.

I use QGIS 3.0.1.

Here's a photo of the attribute table:

And a photo of what it looks like without coloring the different attributes:

Also, here is the actual file:
https://drive.google.com/file/d/1uftFMIc8AVLBWcyV_I5ryn0psw1Y80oM/view?usp=sharing


texts: Set the minimum and maximum values for an elevation raster in R

texts:
I have a raster that displays elevational values at a global scale. I would like to set the minimum and maximum values based on the values I would like to choose.

> elev
class       : RasterLayer
dimensions  : 18000, 43200, 777600000  (nrow, ncol, ncell)
resolution  : 0.008333334, 0.008333334  (x, y)
extent      : -180, 180, -60, 90.00001  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
data source : E:\Data\Elevation_Data\alt\alt
names       : alt
values      : -431, 8233  (min, max)
attributes  :
     ID COUNT
 from: -431     2
 to  : 8233     1

> image(elev, zlim=c(0,1200))

The above command simply plots an image with the necessary values, but I would like a raster with just those limits.

Any thoughts?


texts: Why can QGIS open MrSid files, but GDAL cannot?

texts:
I am able to convert MrSid to GeoTiff in the QGIS GUI, but I would like to do this programmatically through gdal-translate. However GDAL won't recognize the data type, MrSid. But it can when using qgis directly.
gdal_translate -of GTiff -if MrSid ortho_1-1_1n_s_pa079_2015_1.sid /private/var/folders/bk/2m7_mv8d7n5cf43lw0m71j3m0000gn/T/processing_171fdbca8a2e478683a3f68eba3f7141/4d384681ae18421c8d9285b74d514945/OUTPUT.tif
Error Message :
Warning 1: MrSid is not a recognized driver ERROR 4: ortho_1-1_1n_s_pa079_2015_1.sid' not recognized as a supported file format.


texts: CentOS 7, Mapnik, and Tilestache: Seems like image assets are ignored by Mapnik's point/pattern symbolizers

texts:
[Update 6.6.2019]

Ok...new detail. I checked a different Apache error log and discovered that it's littered with these errors, generically corresponding to..

failed to allocate png_ptr

Actual errors (I tried changing the PNG I was using)..

Mapnik LOG> 2019-06-06 12:08:04: Exception caught while loading: '/etc/tilestache/styles/images/click_marker_sm.png' (failed to allocate png_ptr)
libpng warning: Application built with libpng-1.6.16 but running with 1.5.13
Mapnik LOG> 2019-06-06 12:08:04: Exception caught while loading: '/etc/tilestache/styles/images/click_marker_sm.png' (failed to allocate png_ptr)
libpng warning: Application built with libpng-1.6.16 but running with 1.5.13
Mapnik LOG> 2019-06-06 12:08:04: Exception caught while loading: '/etc/tilestache/styles/images/click_marker_sm.png' (failed to allocate png_ptr)

Yikes. Is this a build issue? When Googling it, coincidentally, I find where I asked a largely ignored question over at Unix.SE about this exact Mapnik build. I wonder if the gremlin is hiding somewhere in the weeds behind this comment..

make test ## Mine had a few errors, but I ignored them.. :/

So the error notes Mapnik wants to use libpng-1.6.16, but if I run this..

libpng-config --version
I get.. 1.5.13

Is it maybe enough to just "fix" my default libpng to sync the versions?

[Original Post]

While migrating to a new system (from Ubuntu to CentOS) I've run into an issue where our new Mapnik 3.0.22—built from source—is mostly working, but with one exception. Any layer consuming an external image asset isn't working—specifically, they seem to be rendering 100% alpha/transparent. Any PointSymbolizer or PolygonPatternSymbolizer referencing a PNG file perform the same.

In Tilestache, all of our layers are configured to render from Mapnik XML declarations—so we shouldn't be bumping into Python bindings issues, as 1) that seems to only affect TextSymbolizer, and 2) XML styling is immune to the problem.

Next, I'm sure our data is accessible with the user credentials provided to Mapnik and that the underlying geodata's projection is correct, because tiles will render if I leave the PointSymbolizer parameters empty and force it to use the Mapnik defaults, which seems to be a black dot. For instance here is a tile rendered from the defaults..

<PointSymbolizer allow-overlap="yes">

Meanwhile, the moment I try to point it at a PNG, I just get an empty tile like this (although it's 256x256 of empty pixels)..

<PointSymbolizer allow-overlap="yes" file="/etc/tilestache/styles/images/click_marker.png">

[Update 6.6.2019] FWIW, this is what it should look like (the server in this case is my personal Ubuntu instance with a similar stack, Mapnik 3.1.0, that I sometimes use to experiment or just check my sanity):

The PolygonPatternSymbolizer behaves similarly. For example, here is a tile rendered by our old Ubuntu server using this XML declaration—note the cartographic hashes:

<PolygonPatternSymbolizer file="/etc/tilestache/styles/images/bw-diag-hatch.png"/>

And here's the same tile rendered by the new server using the same XML declaration—note the absence of the cartographic hashes:

I'm not getting any errors pushed up through Apache into the logs. It just seems like Mapnik is ignoring the image assets. Is there something I can do to easily and definitively test whether my Mapnik install will properly consume PNG assets?

I'm to the point where I am questioning my build. :/ Does anyone have any insight or other troubleshooting ideas I can pursue?


texts: QGIS Python composer export problem when using a function

texts:
I try to write a python script in order to do the following:

load a vector layer
open an existing composition
export the composition as PDF

If I use this simplified code in the Python console of QGis.  Everything is working fine, the new layer is added and exported with the composition (you must have at least one composer in your project) :

from PyQt4.QtGui import *
vectorLyr = QgsVectorLayer("d:/test.shp", 'test','ogr')
if not vectorLyr.isValid():
    print "Layer %s did not load" % vectorLyr.name()

vectorLyr.updateExtents()
iface.mapCanvas().refresh()
reg = QgsMapLayerRegistry.instance()
reg.addMapLayer(vectorLyr)
#load first composition
comp = iface.activeComposers()[0].composition()
comp.refreshItems()
comp.setUseAdvancedEffects(True)
comp.exportAsPDF('d:/tmp/test.pdf')
# export is correct

but once I try to generalize the code by creating an export function, the layer is no longer exported with the composition:

from PyQt4.QtGui import *

def export(layerPath, layerName):
    vectorLyr = QgsVectorLayer(layerPath, layerName,"ogr")
    if not vectorLyr.isValid():
        print "Layer %s did not load" % vectorLyr.name()

    vectorLyr.updateExtents()
    iface.mapCanvas().refresh()
    reg = QgsMapLayerRegistry.instance()
    reg.addMapLayer(vectorLyr)
    comp = iface.activeComposers()[0].composition()
    comp.refreshItems()
    comp.setUseAdvancedEffects(True)
    comp.exportAsPDF('d:/tmp/test.pdf')
    return

export("d:/test.shp", "test")
#the vector layer added dynamically is not on the map

I'm new to Python and new to QGis API.  Any idea of what I'm missing?


texts:
I have two shapefiles: one containing a single polyline feature, and the other containing a number of polygons. I'm trying to use ArcPy to identify which polygons are crossed by the line, and then print some of their attributes.

Is 'intersect' the best process for this, or is there a better tool? Everything I'm reading about 'intersect' is suggesting that it will pull points of intersection as opposed to entire polygons.


texts: GDAL python cut geotiff image with geojson file

texts:
I wrote the function which cutting geotiff image by bounding box.

First image is original. At second you can see result off my code. I don't use gdalwarp or any console utilities. But I have no idea how to cut by geojson file. Also I can use only GDAL and numpy modules.

Here is my code:

import os, sys
from osgeo import gdal,gdalconst,osr

def cut_by_bounding_box(min_x, max_x, min_y, max_y):
    xValues = [min_x, max_x]
    yValues = [min_y, max_y]
    #
    # Register Imagine driver and open file
    #
    driver = gdal.GetDriverByName('GTiff')
    driver.Register()
    dataset = gdal.Open(filename)
    if dataset is None:
        print 'Could not open ' + filename
        sys.exit(1)
    #
    # Getting image dimensions
    #
    cols = dataset.RasterXSize
    rows = dataset.RasterYSize
    bands = dataset.RasterCount
    #
    # Getting georeference info
    #
    transform = dataset.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = -transform[5]
    #
    # Computing Point1(i1,j1), Point2(i2,j2)
    #
    i1 = int((xValues[0] - xOrigin) / pixelWidth)
    j1 = int((yOrigin - yValues[0]) / pixelHeight)
    i2 = int((xValues[1] - xOrigin) / pixelWidth)
    j2 = int((yOrigin - yValues[1]) / pixelHeight)

    new_cols = i2 - i1 + 1
    new_rows = j2 - j1 + 1
    #
    # Create  list to store band data in
    #
    band_list = []
    #
    # Read in bands and store all the data in bandList
    #
    for i in range(bands):
        band = dataset.GetRasterBand(i+1) # 1-based index
        data = band.ReadAsArray(i1, j1, new_cols, new_rows)
        band_list.append(data)

    new_x = xOrigin + i1 * pixelWidth
    new_y = yOrigin - j1 * pixelHeight
    new_transform = (new_x, transform[1], transform[2], new_y, transform[4], transform[5])
    #
    # Create gtif file
    #
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(output_file, new_cols, new_rows, 3, gdal.GDT_Byte)
    #
    # Writting output raster
    #
    for j in range(bands):
        data = band_list[j]
        dst_ds.GetRasterBand(j+1).WriteArray(data)
    #
    # Setting extension of output raster
    #
    dst_ds.SetGeoTransform(new_transform)
    wkt = dataset.GetProjection()
    #
    # Setting spatial reference of output raster
    #
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    dst_ds.SetProjection( srs.ExportToWkt() )
    #
    # Close output raster dataset
    #
    dataset = None
    dst_ds = None

if __name__ == '__main__':
    # Imput/output file name and set directory
    os.chdir('/home/sant/test/satellite_images')
    filename = '20160501.tif'
    output_file = '/home/sant/test/20160501_cutted_by_bounding_box.tif'
    cut_by_bounding_box(531961.73, 535987.34, 4894164.57, 4888631.61)
    print 'cutter.py script done!'


texts: How to change default symbol in Categorized symbol renderer (python console) or use different icons as vector layer symbols?

texts:
I have a point shapefile in QGIS 2.10 and I want to display one of its fields attributes based on its unique values (I always know that it has 4 unique values and their values so I do not care for the moment to extract them) I just want points to be displayed with different symbols depending on its value for the specified field.

The categorized symbol renderer styles features based on unique values in the layer's attribute table (QGIS python programming cookbook, chapter 5, page: 145). Looking at the API for QgsSymbolV2 I think that maybe I have to use the public member fuction changeSymbolLayer  instead of

QgsSymbolV2.defaultSymbol(lyr.geometryType())

but I get very confused on how to follow the API.

Then, there is the option to use icons as vector layer symbols (QGIS python programming cookbook, chapter 5, page: 140) but I am not sure in how to achieve to have different icons for different categories.


texts: Convert Date Field to String Field Arcpy

texts:
I have a Date field in an ArcGIS file geodatabase feature class called "DateTemp1." I created a new text field called "DateTemp2" and want to calculate that based on the DateTemp1 field, but I am getting a unicode error.

TypeError: descriptor 'strftime' requires a 'datetime.date' object but received a 'unicode'

 arcpy.CalculateField_management("all_reports.gdb/lsr_reports1", "DateTemp2", "datetime.datetime.strftime(!DateTemp1!, '%Y%m%d')", "PYTHON_9.3", "")


texts:
I have some problem with my script.

I want to add my picture link. I save the picture link at some table in geoserver

this is my script


texts: How to transform Lambert coordinates in C#?

texts:
I want to transform Lambert (EPSG:3035) coordinates to EPSG:4326 in C#.

I have tried this:

var lambert = cf.CreateFromWkt("PROJCS[\"ETRS89 / ETRS-LAEA\",GEOGCS[\"ETRS89
\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS
1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6258
\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree
\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4258
\"]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center
\",52],PARAMETER[\"longitude_of_center\",10],PARAMETER[\"false_easting
\",4321000],PARAMETER[\"false_northing\",3210000],UNIT[\"metre\",1,AUTHORITY[\"EPSG
\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"3035\"]]");

Error: Projection Lambert_Azimuthal_Equal_Area is not supported.

As far as I understand, ProjNet should support Lambert. Can someone help me? Or show me another way how to achieve my goal?


texts: Which projection casts rays from a point on the surface of the globe to a cylinder?

texts:
What is the name of a map projection that casts rays from a point on the surface of the globe to a cylinder?

Some people who make sundials put a map of Earth on the sundial.
I want to use the "right" map projection such that, when the nodus casts a shadow on a particular city on the map, the sun is directly overhead that city on the real Earth.

If the sundial is in the form of (part of) a cylinder, and the nodus fixed to a point on the axis of the cylinder, then such a map uses the cylindric perspective projection.

If the sundial is entirely composed of flat surfaces, then such a map uses a gnomonic projection on each surface.

A few sundials are in the form of (part of) a cylinder, and the nodus is a small notch or hole in the curved surface of the cylinder. Carl Sabanski calls this a "Cylinder Sundial: Type 1".
Is there a name for the kind of map projection I want to use in this sundial?

In the form of an analogy question:
gnomonic projection : azimuthal conformal projection :: cylindric perspective projection : __?

Edit:
Oops.
My description seems to be a mishmash of 2 different projections.
For lack of better names, for now I'm going to call them the "point-on-globe-in-cylinder projection" and the "Type 1 cylinder sundial projection".
I hope that I will soon learn some other, less awkward name for each projection.

After posting the original question, I realized that the "point-on-globe-in-cylinder projection" is not the one I want to put on my "Type 1 cylinder sundial".

The "point-on-globe-in-cylinder projection" involves a globe inscribed in a cylinder, where both of them have radius R, like most cylindrical projections.
Like the azimuthal conformal projection, there is very low distortion around one special center point on the map.
Like the azimuthal conformal projection, rays start at the antipode of that center point, then hit some city on that globe, travel a bit further and hit the paper where that city is mapped.
Like most cylindrical projections, since the globe touches a circumference of the cylinder along some great circle, the map is pretty accurate at points near that great circle.
With a very long cylinder, almost all the globe can be mapped.

The "Type 1 cylinder sundial projection" is the projection I want to place on my "Type 1 cylinder sundial".
The globe has radius R, while the cylinder has radius R/2.
The globe is positioned so the geometric center of the earth is at the nodus hole on the surface of the cylinder.
The point on the surface of the cylinder diametrically opposite the nodus is a special central point.
The surface of the cylinder is tangent to the surface of the globe (only) at that central point.
Along the straight line from that central point along the length of the cylinder, it's approximately a gnomonic projection -- rays start from the center of the globe, then hit a city on the globe, then travel a bit further and hit the paper where that city is mapped.
Along the circumference from that central point around the cylinder to the nodus hole, rays start from the center of the globe, then hit the paper cylinder where a city is mapped, and then hit that city on the globe.
With a very long cylinder, almost half the globe can be mapped, the hemisphere around the central point.
Since there is a small hole punched in the surface of the cylinder to form the nodus, rays from the center of the nodus to the antipode of the central point (and the hemisphere around that antipode) never strike the cylinder, and so are not mapped anywhere on the cylinder.

Is there a better name for what I'm calling the "point-on-globe-in-cylinder projection"?
Is there a better name for what I'm calling the "Type 1 cylinder sundial projection"?

Should I split this question up into 2 different questions,
each one asking about 1 of these 2 projections?


texts:
Using Arcpy, how can I determine min and max elevation from a raster DEM layer within my current extent? I'm currently using Data Driven Pages to centre the layout on each feature extent.

I need to determine this elevation difference value in order to get Python to automatically activate the appropriate contour interval feature class (25m, 100m, etc)

The current logic I'm thinking of is to select all cells within current extent, (what is the most convenient way, or arcpy methods() to do this?) place their values in a list, then:

minElev = 99999
maxElev = -99999

for i in myList:
    if i < minElev:
        minElev = i

for i in myList:
    if i > maxElev:
        maxElev = i

then subtract max and min to get the difference value.

I don't think this is the easiest way, is there someone experienced who knows a faster way?


texts: What USGS maps exist in the 7.5' × 15' metric 1:25,000 series?

texts: While browsing through the library, I discovered that USGS has, in fact, produced some metric maps at 1:25,000.  One example is Thendara (NY), 43075-F1-TM-025.  As far as I can tell, those maps are sporadic.  I do find some at the USGS website.  What is the nature of this series?  Is it a past project since abandoned, is it only for some areas, and does it exist in place of or in parallel to the imperial 1:24,000 series?  What areas are covered by those maps?


texts:
When I registered a database data store to ArcGIS Server in ArcCatalog, the following error occurs:

  Underlying DBMS error[ORA-12154: TNS:could not resolve the connect
  identifier specified No extended error.]

I set up the ArcGIS Server, ArcGIS for Desktop, and Oracle are all installed on one machine.
ArcGIS for Server 10.2.2 64bit
ArcCatalog 10.2.2 32bit
Oracle Instant Client 11.2.0 64bit
Oracle Database 11g Enterprise Edition 11.2.0 32bit

I have successfully made the Oracle 11g "Enabled Enterprise Geodatabase" through the ArcCatalog, and it seems nothing wrong.

Database Platfrom: Oracle
Instance: orcl
Authentication Type: Database authentication
User name: gesbar2
password: gesbar2

And I have added the Oracle Instant Client directory to the PATH and environment variable, I am sure that the 64-bit DLL is found before the 32-bit one.

But i still get the error message.


texts: Repoint feature layers from one SDE to another

texts:
I am trying to repoint layer (.lyr) files to a new SDE connection. In a test project folder I have two test layer files.

When using lyr.findAndReplaceWorkspacePath() the script runs with no errors, but only the first layer in the directory gets properly updated to point to the new SDE path.

When using lyr.replaceDataSource()the script produces an error while attempting to repoint the second layer to the new SDE path:

ValueError: Layer: Unexpected error

(Error is referring to the line containing lyr.replaceDataSource())

So something is running afoul on the second iteration of the for loop.

I have tested multiple .lyr files, and it's nothing wrong with any of the actual .lyr files. It's just whichever .lyr file is up second. So whichever is first, is properly repointed to SDE, and whichever is second, is not.

import arcpy, os

layersFolder = r"I:\FolderX\Kristen\FolderY\Repoint_lyrs\TestLayers\MultipleLayers"
newSDEpath = r"I:\FolderX\Kristen\FolderY\Repoint_lyrs\gis_sde.sde"

for dirname, dirnames, filenames in os.walk(layersFolder):

    for filename in filenames:
        path = os.path.join(dirname, filename)
        basename, extension = os.path.splitext(path)

        if extension == ".lyr":
            lyr = arcpy.mapping.Layer(path)

            if lyr.isFeatureLayer:

                if lyr.supports("DATASOURCE"):

                    lyrName = lyr.datasetName
                    print lyrName

                    print '\t' + "newSDEpath: " + newSDEpath
                    print '\t' + "lyr.datasetName: " + lyr.datasetName

                    lyr.replaceDataSource(newSDEpath, "SDE_WORKSPACE", lyr.datasetName, True)

                    #lyr.findAndReplaceWorkspacePath(lyr.dataSource, newSDEpath, True)

                    lyr.save()

                    print '\t' + "Changed path: " + lyr.dataSource + '\n'

                    del lyr


texts: Creating new layer using selection from attribute table in QGIS?

texts:
How do I create a new layer using a selection from an attribute table in QGIS 2.10.1?

Previous answers on this topic do not help -perhaps because they relate to older versions of QGIS


texts: Clip the contour map with the polygon using R software?

texts:
I am wondering how can I create the contour lines from the discrete x, y and z dataset and then add the shapefile to cut the contour lines by polygon.

I want to clip the contour map by the shapefile. I have created the contour map and added the shapefile but I am not able to cut the shapefile with polygon.

I used akima package to create the contour map. I have posted this question in stackoverflow R help but I didn't get any help so far and I desperately need help on this issue.

I will post my code and data links here:

The script file image.scale.R can be found on the following location "https://www.dropbox.com/s/2f5s7cc02fpozk7/image.scale.R "

The code that I have used so far is as follows:

## Required packages
library(maptools)
library(rgdal)
library(sp)
library(maptools)
library(sm)
require(akima)
require(spplot)
library(raster)
library(rgeos)

## Set Working Directory
setwd("C:\\Users\\jdbaba\\Documents\\R working folder\\shape")

## Read Data from a file
age2100 <- read.table("temp.csv",header=TRUE,sep=",")

x <- age2100$x
y <- age2100$y
z <- age2100$z

####################################
##Load the shapefile
#####################################
shapefile <- readShapePoly("Export_Output_4.shp")

fld <- interp(x,y,z)

par(mar=c(5,5,1,1)) filled.contour(fld)

###Import the image.scale
source source("image.scale.R")
# http://menugget.blogspot.de/2011/08/adding-scale-to-image-plot.html

x11(width=8, height=7)
layout(matrix(c(1,2), nrow=1, ncol=2), widths=c(6,1), height=6, respect=TRUE)
layout.show(2)

par(mar=c(4,4,1,2))
image(fld,axes=T)
contour(fld, add=TRUE)
#points(age2100$x,age2100$y, pch=".", cex=2,legend=F)
plot(shapefile,add=T,lwd=2)
box()
par(mar=c(4,0,1,4))
image.scale(fld, xlab="Eastings", ylab="Northings", xaxt="n", yaxt="n", horiz=FALSE)

axis(4)
mtext("Salinity", side=4, line=2.5)

Now, I want to get rid of the colored gradients and the contours from the polygon shapefile and only leave the intersection part.

Any help is highly appreciated.


texts:
I have been using GDAL/OGR in a project and it is going quite well. I have been re-projecting data from various source formats into a single CRS (using the C# bindings), but one issue I have come across is with the associated OGR style information.

Most of the font sizes in the source documents are stored in "map ground units", so if I just copy the OGR style directly the text size will be scaled up or down depending on the differences in the map units.

So my question is, what is the best way of converting the map ground units in the style specification from one CRS to another, or in pixel units?

Many thanks for any help you can provide,


texts:
Is it possible to have the legend scroll using the OpenGeo Suite SDK and GXP. Following this tutorial Tutorial LinkIt is possible to have it any size if external as the height and width can be changed but when embedded into the west container, if too many layers are selected the legend doesn't fully display each layer?

Suggestions or ideas?


texts: ArcGIS Pro freezes when using Measure tool

texts:
Every time I finish measuring the length of a line feature (specifically 500 m sections of a river) using the Measure tool in the Map ribbon, ArcGIS Pro freezes up and the swirling blue lagging cursor appears next to my mouse cursor. It takes several minutes for it to recover and in the meantime, I cannot do anything whilst I wait, including using the 'Clear Results' eraser button.
The line features have a spatial index 'FDO_Shape'
Attribute INdex: 'FDO_OBJECTID'
Data type: File Geodatabase Feature Class
Feature Type: Simple
Geometry Type: Line
This is really affecting my workoutput and I think there must be a bug in the system affecting the Measure tool.
Can anyone help fix this issue so that I can use the Measure tool without it freezing up my project?


texts: Clip raster (DEM) with vector layor in QGIS

texts:
I am trying to trim/clip a DEM with a shapefile layer. I read somewhere that the GDAL should be able to get the job done, but can't seem to get it to work.

I have a raster layer and a shapefile layer, then using the "clip raster by extent" where the extent is generated from the shapefile layer.

The red square is the cutting line, and the black and white background is the DEM maps which I want to cut/trim

This is the input into the dialogue

The elevation data needs to be kept as the map will later be used to create a 3D model of the area.


texts: How to join attributes of two different road networks with different segment lengths

texts:
I ran an accessibility analysis on a road network with GRASS GIS (v.net.iso Tool) from different nodes in the network.

Using QGIS 2.18, I need to attach the resulting layer to the road network again. This is because v.net.iso puts out a layer containing just the distance categories. In my case 3 categories since I ran an analysis with 2 different distances (5km, 7,5km and above is cat 3).

Now if I use the "Join attribute by location" tool, I can choose the starting layer to which I want to join the attributes of the road network layer (street numbers, car counts etc.). In my case this would be the layer which GRASS put out.

The problem is, that with "Join attribute by location" a part of road B, which belongs to distance category 1 is attributed with features of road A because this is, where category 1 "started". I tried to show this in the picture below. I would need to "cut" the roads layers with its segments with the category layer but I dont't know how.


texts:
Is there a way in SciPy or NumPy to buffer values in an array?

I have several rasters that I read as arrays using gdal to do some math/masking and then I write the final array back to geoTIFF. One of the layers I use to make my outputs needs to have the 0's in the array be buffered by one pixel before it can be used as a mask. Is there a way to do this to an array in Python? So if my input array looked like this:

1 0 1 1 1 1
0 0 1 1 1 1
1 0 1 1 1 1
1 1 1 0 1 1
1 1 1 1 1 1

Then the out array would look like this:

0 0 0 1 1 1
0 0 0 1 1 1
0 0 0 0 0 1
0 0 0 0 0 1
1 1 0 0 0 1

If there's not a straightforward way to do this with arrays, what's the best way to do this in GDAL? I only want the 0's in the raster/array to be buffered with 0's, and just by one pixel.


texts:
I wanted to create a new SRID for use in PostGIS. My goal is to have a SRS that can cover an entire state (instead of have multiple zones as in SPCS). For example, I wanted to create a new SRS for Kansas, based on a compromise between

Kansas North (EPSG 3427):

  PROJCS["NAD83(HARN) / Kansas North (ft
  US)",GEOGCS["NAD83(HARN)",DATUM["NAD83_High_Accuracy_Regional_Network",SPHEROID["GRS
  1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6152"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4152"]],UNIT["US
  survey
  foot",0.3048006096012192,AUTHORITY["EPSG","9003"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",39.78333333333333],PARAMETER["standard_parallel_2",38.71666666666667],PARAMETER["latitude_of_origin",38.33333333333334],PARAMETER["central_meridian",-98],PARAMETER["false_easting",1312333.3333],PARAMETER["false_northing",0],AUTHORITY["EPSG","3427"],AXIS["X",EAST],AXIS["Y",NORTH]]

Kansas South (EPSG 3428):

  PROJCS["NAD83(HARN) / Kansas South (ft
  US)",GEOGCS["NAD83(HARN)",DATUM["NAD83_High_Accuracy_Regional_Network",SPHEROID["GRS
  1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6152"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4152"]],UNIT["US
  survey
  foot",0.3048006096012192,AUTHORITY["EPSG","9003"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",38.56666666666667],PARAMETER["standard_parallel_2",37.26666666666667],PARAMETER["latitude_of_origin",36.66666666666666],PARAMETER["central_meridian",-98.5],PARAMETER["false_easting",1312333.3333],PARAMETER["false_northing",1312333.3333],AUTHORITY["EPSG","3428"],AXIS["X",EAST],AXIS["Y",NORTH]]

I wanted to combine these into one coordinate reference system. My questions are:

How can I create a new SRS based on OGCWKT?
(There's a related question Create custom coordinate system in postgis, but it turns out that the OP is just wasn't able to find an existing SRID, and the answer is about finding and inserting an existing SRID into PostGIS).

My plan is to use standard_parallel_1 from the North Zone, and standard_parallel_2 from the South Zone to create the new SRS.

Can I just make a copy of the North Zone, and change the two standard parallels to get the new OGCWKT, or are there additional parameters I must change?

I have basic knowledge about projections in general, but not a lot.


texts:
I would like to apply Moran's I spatial auto-correlation measure to 2D arrays or raster files in Python and I am wondering if there is an available implementation for it.

I know Pysal offers an implementation but I was not able to figure how to apply it to 2D data.

Any other libraries/implementations that would allow me to do that?


texts: DXF (points and lines) in QGIS - Doesn't import z coordinates

texts:
Trying to get to grips with QGIS, coming from a MapInfo and Civil3D background.

I have a DXF file with points (X,Y,Z) and lines.  I want to open this in QGIS and use the interpolate plugin to create my DEM.

When I open the DXF I assign a coordinate system, however the Z attributes don't seem to come across and so the interpolate function fails.

I could be doing something wrong on a very basic level - any pointers on where to start?


texts: What's the name for small area statistics in France?

texts: US has Census blocks and block groups. Canada has dissemination blocks and dissemination areas. UK has output areas and super output areas. Is there an equivalent in France—statistical areas defined as having, at most, several thousand people?


texts:
Having trouble upgrading from QGIS 1.8, which ran just fine, to QGIS 2.8.

I'm running 64-bit Windows 7. After uninstalling QGIS 1.8, I've been trying for several days to install the 32-bit version of the QGIS 2.8 Standalone Installer.

The installer runs, but there's no program to start. The Start Menu folder is empty, and there is no qgis.bat in the directory that I installed to (C:\qgis). Trying to run qgis-bin.exe throws an error telling me that qgis_core.dll can't be found. The postinstall.log file has 15 instances of the error "System cannot find the path specified."

Any ideas?


texts: Making polygons out of X,Y-point groups?

texts:
I have (X,Y)-coordinates for different places in a single file. How can I make polygons from a group of points separated by NULL? A common method 'Points to line' does not seem to work.

For instance, from the picture above the first six points ("SN" from 1 to 6) will form the first polygon, next four points ("SN" from 8 to 11) will form the second polygon and so on.


texts: Converting HDF file (MODIS product) in ArcGIS as correctly georeferenced raster by ENVI?

texts:
Which algorithm I should use in ENVI 5.2 to open the MODIS product (MOD13Q1 or MODIS/Terra Vegetation Indices 16-Day L3 Global 250m Grid SIN V006), given as  HDF file to ArcGIS 10.5 as correctly  georeferenced raster showing the value of NDVI index?

The product is described here:
https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod13q1_v006

I tried to do it by  ENVI 5.2. I opened it there (Open as - EO-1 - HDF4), but  when I exported it directly by operation “Send to ArcMap”, it was opened in ArcMap with incorrect projection. The same things happened when I imported this hdf file to ASCII and after that in ArcMap used the “ArcCII to Raster” Tool. The resulting raster has incorrect projection.
In the description of this product (by link above) says that it has Sinusoidal Coordinate system. So I suppose I should change it in the properties of file. However, I can’t change this coordinate system in the files, which I created by ENVI.
So, I need to know – which algorithm I should use in ENVI 5.2?


texts: How to use ST_Union with polygon and multipolygon geometry types?

texts:
I am attempting to combine the geometries for multiple records using the following:
INSERT INTO schema.boundary (name, geom)
SELECT 'Name' ,  ST_Union(geom) FROM schema.boundary WHERE id = 0;

This returns the error:
ERROR:  Geometry type (Polygon) does not match column type (MultiPolygon)
********** Error **********

So I then tried the following:
INSERT INTO schema.boundary (name, geom)
SELECT 'Name' ,  ST_Union(ST_Multi(geom)) FROM schema.boundary WHERE id = 0;

This returns the same error. I managed to combine all of the target geometry fields except for one. Checking the geometry type using ST_GeometryType(geom) for the all records returns 'ST MultiPolygon'. I've also tried updating the single record using ST_Multi(geom).
Is there another PostGIS function I can use to combine incompatible geometry types?


texts: How do you pass parameters for Proj4js implementation with OpenLayers

texts:
I am using Proj4js along with OpenLayers for the transformation of x,y (projection type is Transverse Mercator) to long, lat.

While executing "inverse()" method what I noticed that in tmerc.js some "this.ep2" parameter is used with in the formula and the value of this parameter is "NaN".

Because of this the transformation is not working.

I am guessing that this is happening because perhaps I am not initializing required parameter with in the definition file of the source coordinate system.

So, could anyone help me in knowing which parameter for "this.ep2" has to be passed in definition file?


texts: Transparency setting for polygon borders in ArcMap

texts:
I'm working on a map using ArcGIS Desktop 10.8 and need to create a style with 50% borders transparency which I can then apply to different drawing elements like polygons, lines and text. Similar to the yellow elements in this BBC map.

How can I achieve this?


texts:
I'm looking for very high resolution stereo imagery to perform a DEM extraction. The only place I can find sample data is http://www.digitalglobe.com/product-samples but the form does not work for me.

Any idea where I can find the data I'm looking for? QuickBird 2, IKONOS, WorldView, GeoEye, RapidEye, or some other very high resolution source and of mountainous terrain would be ideal.


texts:
I am looking for an easy way to convert a polygon to a polyline in r.
I have tried this but dit not succeed:

polygon <- readOGR("~path/", "polygon.shp")
polyline <- as(polygon, "SpatialLinesDataFrame")
writeOGR(polyline, dsn="~path/", layer="polyline.shp", driver="ESRI Shapefile", overwrite_layer=TRUE)


texts: What causes ArcPy RuntimeError: TextElementObject: Error in setting text?

texts:
I am trying to generate a map book with the following code that I downloaded but I keep getting an errror:

  *File "Y:\Planning\GIS Planning Analyst\projects\county_map_book\9_GenerateMapBook.py", line 18, in
  <"module>
      titleText.text = title   File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\arcobjects_base.py", line 77, in
  _set
      return setattr(self._arc_object, attr_name, ao) RuntimeError: TextElementObject: Error in setting text*

#Read values from input dialog
title = arcpy.GetParameterAsText(0)
organization = arcpy.GetParameterAsText(1)

arcpy.AddMessage("    Updating Title Page ...")
titleMXD = arcpy.mapping.MapDocument("path\county_map_book\MapBook_TitlePage.mxd")
titleText = arcpy.mapping.ListLayoutElements(titleMXD, "TEXT_ELEMENT", "Title")[0]
titleText.text = title
arcpy.mapping.ExportToPDF(titleMXD,"path\county_map_book\Output\1_TitlePage.pdf")
del titleMXD


texts: Seeking commercial web service for LIDAR/DSM data in USA

texts:
I apologize for not knowing all the proper terminology to ask this question. I am seeking a commercial service for providing pay-per-use high resolution elevation data in North America.  It must come with defined ToS (terms of service -- a guarantee of availability).  A paid service is best.

The ideal service would return a chunk of elevation data, perhaps in ARC/INFO ASCII GRID format, on demand.  Ideally this grid would be seamlessly centered on a given coordinate.  For this application I need a regular grid of elevation points perhaps 2000 or 3000 grid elements in the x and y direction (if the data returned were richer LIDAR data, I would need to interpolate into a grid using las2dem or some such).

I am aware of scattered sources of free LIDAR elevation data:

http://lidar.cr.usgs.gov/
http://www.opentopography.org/
http://calm.geo.berkeley.edu/ncalm/ddc.html

I'm seeking something that is (a) offering wider coverage with less hassle (b) a web service returning data in close to real time (c) supplying seamless data centered on a point on earth, rather than tiles (d) commercially usable and supported with terms of service (e) command line scriptable (e.g. using wget or other scriptable tools outside of a GIS).  For example:

wget http://ppv.example.org?lat=37.1234&lon=-122.1234
&rows=1024&columns=1024&format=ascii_esri_grid&data_type=dsm:best_available

Is this a crazy desire, or is there something that fits the bill?

Possible question keywords: WFS, web service,hosted LIDAR data server


texts: How to extract values from a raster stack within overlapping buffers?

texts:
In R, I have a raster stack with 6 variables. I would like to extract the mean and maximum value of each variable at 300 coordinates within a 500m buffer radius and a 1km buffer radius.

In ArcGIS 10.2 there are problems extracting variables from buffers if there is any overlap. Is this a problem in R?

If not, how do you make buffers around a point in R?
And how would I extract the mean and maximum value at the two scales surrounding each point?


texts: Adding HTML to folium popup in a for loop

texts:
I want to put HTML formatting into a folium map popup. How do I format the html popup inside my for loop?

My code:

import folium
import pandas
from folium import IFrame

data=pandas.read_csv("addresses.txt")
lat=list(data["LAT"])
lon=list(data["LON"])
name=list(data["NAME"])
desc=list(data["DESCRIPTION"])

map=folium.Map(location=[40.748817,73.985429],zoom_start=14,tiles="stamenterrain")
#folium.TileLayer('cartodbpositron').add_to(map)

for lt,ln,nm,ds in zip(lat,lon,name,desc):
    test = folium.Html('nm
desc', script=True)
    fg.add_child(folium.CircleMarker(location=[lt,ln],radius=6,color='grey',fill_color='yellow',
popup=test))

map.add_child(fg)
map.save("Map1.html")


texts: Assigning a variable for 'iface.activeLayer()' when starting QGIS

texts:
I use iface.activeLayer() intensively. Therefore, when I open QGIS, firstly I use the following script in QGIS Python Editor not to enter iface.activeLayer() everytime.

from qgis.utils import iface
L = iface.activeLayer()

def activeLayerChanged():
    L = iface.activeLayer()

iface.currentLayerChanged.connect(activeLayerChanged)

I thought it would work if I added it to startup.py file. I tried to add the script to startup.py but I got error: name 'L' is not defined.

How can I assing iface.activeLayer() to L (or any single-letter variable)  when starting QGIS, so that everytime I enter L, I can get the active layer?


texts:
New to JavaScript, coming from Python background. Created a custom control and have added buttons inside the _div using a for loop and innerHTML.

I want to zoom to and open the popup of a layer of the specific button I click. When I click on the button the function zooms to the marker, but does not open the popup.

The popup is definitely attached, since when I click on the marker the popup opens. When I enter the same command into the web console tool, the map zooms to and then opens the popup. Not sure what I am doing wrong.

What really confuses me is the fact that the function behaves differently when I use the button vs when I use the web console.

//variables
var map_object = L.map('mapid');
var map_tile_object = L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png');
var resume_markers = []
var resume_control = L.control();
var resume_control_title = '
Interactive Map
';
var resume_control_description = '
Click on item for more information
';

//functions
function zoomTo(num) {
  var feature_cord = resumeData[num].geometry.coordinates;
  var markerBounds = L.latLng(feature_cord[1], feature_cord[0]);
  map_object.setView(markerBounds, 18);
  console.log(num);
  resume_markers[num].openPopup();
};

function main(map, tiles, data, resume, title, description){
    map.setView([41.505, -77], 6);
    tiles.addTo(map);
    for (x in data){
      var popup_content = data[x].properties.responsibilities;
      var marker = L.geoJSON(data[x]).bindPopup(popup_content);
      resume_markers.push(marker);
    };

    for (x in resume_markers){
      resume_markers[x].addTo(map);
    };

    resume.onAdd = function(map) {
      this._div = L.DomUtil.create('div', 'resume');
      this._div.innerHTML = title + '
' + description;
      return this._div;
    };

    resume.addTo(map);
    for (var i = 0; i < resume_markers.length; ++i){
      resume._div.innerHTML +=
      '<button onclick=zoomTo('+ i + ')>' +
      data[i].properties.position +
      '</button>';
    };

};

//execute
var app = main(
  map_object,
  map_tile_object,
  resumeData,
  resume_control,
  resume_control_title,
  resume_control_description
);


texts: Importing .shp Files Into R: readOGR error

texts:
I'm trying to import shapefiles into R but receiving errors. I'm using rgdal package and doing the following:

x <- readOGR(".", "SNP Basefile")

  Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding,
  use_iconv = use_iconv,  :  Cannot open layer

I used dsn="." since I had placed the shapefile in the current working directory. "SNP Basefile.shp" is the shapefile I'm trying to import.


texts: Mapserver Mount Error: Failed to resolve server ip.address.of.nas: Name or service not known

texts:
I followed this guide to install mapserver on ubuntu:
https://goworldwind.org/mapserver-and-data-installation/
, and when I set the step 3: Mount the MapServer Configuration and Data, I did

sudo mount ip.address.of.nas:name.of.vbox.share /data

There was an error:

Failed to resolve server ip.address.of.nas: Name or service not known

I used the VMware Workstation, anyone know how to solve this problem?


texts:
I'm fairly new to QGIS but thought I kinda got the CRS issues. I'm using version 2.8.3 on Windows.

However I'm currently trying to import an asc-file from the German weather service.  The problem is that the layer has- once imported- a terrible offset.  The map with wind-data of Germany is suddenly located over the black sea.

I've tried all things I could find. Predefining the CRS-System, with ->setting->options, assign projection, tried out several different crs, but it all won't work.

The file description from the weather service says it is

  Gauss-Kruger-Projection with reference meridian 9 degree east (Bessel-Ellipsoid)

But I have no idea what to do with this information since I can't find this CRS among the offerer CRS.

Here you can find the files, which are giving me problems. It's the last asc-file in the folder: ftp-cdc.dwd.de


texts:
I'm looking for some free geocoding service with possibility to retrieve contact information to specific place like phone number, email etc.

Currently I'm using MapQuest:
http://www.mapquestapi.com/geocoding/
and they are returning only localization and name


texts: QGIS change pixel size after multiplication in raster calculator

texts:
QGIS change raster size after multiplication in raster calculator:

I have layer with size 25x25 pixel as you can see in the picture:

When I multiple this layer * 0 in raster calculator (I need clean layer) I get layer with 100x100 pixel size:

What's wrong? is it bug?


texts: Calculating absolute precision confidence number from Dilution of Precision indicators from GPS receiver

texts:
End-user GPS receivers typically show an absolute confidence indication, such as Within 25 Metres, for the calculated position.

I am parsing NMEA data from a GPS unit and wan't to calculate a similar quality of service indication. The NMEA data exposes a Dilution of Precision indication (actually, one for both the horizontal and vertical position, HDOP/VDOP).

How do I go from DOP value(s) to a concrete confidence measurement (say, metres)?


texts: Convert contours to polygons maintaining attributes using Python?

texts:
I have a series of polylines(contours) that i want to convert to polygons using Python (GeoPandas, Fiona, Shapely). Because these are contours, there will be donut polygons, can anyone point me to how to do this? Each contour is a polygon, with there is no overlap of polygonsdonut polygons are polygons with holes in them. the holes are other polygons that have been created by inner contours

my non working code

def polyline_to_polygon(input_file, output_file):
    polylines = fiona.open(input_file)
    geom = [x['geometry'] for x in polylines]
    coords = [x['coordinates']for x in geom[:-1]]
    # 'ndviChange': 'float'}} should be carried over from the polyline
    schema =  {'geometry': 'Polygon','properties': {'id': 'int'}}
    with fiona.open(output_file, 'w', 'ESRI Shapefile', schema) as c:
        for poly_id,polygon in enumerate(polygonize(coords)):
            c.write({
                'geometry': mapping(polygon),
                'properties': {'id': poly_id},
            })

This works to generate the correct polygons, but doesn't include the ndvi value in the output shapefile. So there should be a ndviChange value in the polygon. it is in the polyline fc but doesn't carry over.

Hicham, i tried that and received

Must have equal len keys and value when setting with an iterable

gdf.loc[index, 'geometry'] = gdf.loc[index, 'geometry'].difference(gdf.loc[index_, 'geometry'])


texts: Not Able to Add Breaks For Feature Service Graphics

texts:
Using ArcGIS JavaScript API 3.24, can you please take a look at this snippet and let me know why I am not able to add Symbology to Feature Service out put by using ClassBreaksRenderer

Here is the code

<!DOCTYPE html>

<meta charset="utf-8">
<meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no">
<title>FeatureLayer</title>

<link rel="stylesheet" href="https://js.arcgis.com/3.24/esri/css/esri.css">


texts: GeoTiff specs for coordinate transformation definition

texts:
I would like to define a GeoTiff coordinate transformation by a set of GCPs only. The intention is to publish the GeoTiff as a coverage in GeoServer.

However if I specify tie points only, GeoServer complains that I also need to specify the scale tag, which then results in only one tie point used, eg no inclination or shearing. Is this normal? Can't I use GCPs alone for projection?

On the other hand I could define a transform, but the transformation is rather limited compared to let say a set of 10 GCPs.

Thanks for any help,

Thomas


texts: Understanding why Split Polygon tool is greyed out?

texts:
I have a wetlands dataset (blue in the pic below) and a parcel dataset (green). I would like to clip the parcel dataset where it intersects the wetlands dataset, so that the parcel layer will be empty where the wetlands go through.

The Advanced Editing tool 'Split Polygons' seems to be my best option, but it is greyed out. I have an advanced license, so I know that's not the issue.

Does anyone know what could be causing this to happen?


texts:
I am using qgis 2.14.17 with ubuntu16. Having as input data a DEM of the region and a river network.
In the manual of qgis they suggest to use a toll of qgis: Catchment area (Parallel); I dont have this in my toolbox, anybody made the same experience and could tell me the reason for it and how to get it. It seems to offer the best solution for me.
Tools which are available are from GRASS7 and TauDEM, the latter need much more data it seems. And I am not sure if I am able to create these just out of a DEM? The GRASS7 has one tool r.water.outlet: but here one has to put the outlet coordinates, which clearly I could get out of the river network, but the idea is that they are generated by delineating the watersheds.

Somebody has an idea how to solve it?

I managed now to get SAGA run in the toolbox, nevertheless would be great if somebody could help me with the qgis tool, since it should appear?!