texts: Confused about "during daylight" VS "daily" solar radiation

texts:
Since light is important for algal growth, its effect is simulated in many water quality models, e.g.,  WASP model (see Equation 6 in its manual for multiple algae). However, below is  what I cannot understand well about this Equation 6:

Why Ia = 0.9*I/f (where Ia is the average incident light intensity during daylight hours; I is incident solar radiation; f is fraction of day that is daylight)? What does this 0.9 come from?


texts: How to refresh QGIS Relation Reference widget without closing and reopening form?

texts:
The Relation Reference widget has buttons to open a 2nd form so that records in the referenced (parent) layer can be edited or new records added.

However, when the edits are made to the parent layer and 'OK' is clicked on the 2nd form - the edits are not reflected in the list of values in the widget, unless both these are performed:

the parent table is saved
the form containing the Relation Reference widget is closed and reopened

How else can the Relation Reference widget be refreshed ?


texts: Increasing PostGIS/PostgreSQL memory to avoid Out Of Memory (OOM) error?

texts:
This is a follow-up question to my earlier question about the memory management of some PostGIS/PostgreSQL functions:

How do PostGIS functions handle buffers and external storage?

As suggested by the accepted answer, some of the PostGIS functions (e.g. that rely on GEOS or those in pgrouting)

  do not have any special handling for data sets that are larger than
  can fit in memory. If you run them on such data sets, you'll just OOM
  the backend.

So, my question here is that

If I have a machine with lots of memory, how do I config PostgreSQL/PostGIS so that I don't get OOM error for larger datasets?

I see here that there are a few parameters:

work_mem, shared_buffers, maintenance_work_mem

For example, if I have a pgrouting memory taking 100M worth of memory. Which ones in above should be set and what values should I use to ensure I don't get an OOM? (assuming no other concurrent use of the database).


texts: Extruding a 3D stack of polygons on top of a map, e.g. with KML

texts:
I have data for the composition of rock layers in several polygons on earth. I would like to plot a stacked chart of the rock layers, colored based on layer metadata. Here's a mockup, where rock layers are extruded downward from surface polygons:

So far, I've been experimenting with making this visualization using KML. From the samples, I've found how to extrude above a single polygon, but I haven't been able to make a stack.

How do I modify this KML to produce a stacked set of extruded polygons? The two extrusions are unfortunately intersecting currently.

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Document>
    <name>KmlFile</name>
    <Style id="transBluePoly">
      <LineStyle>
        <width>1.5</width>
      </LineStyle>
      <PolyStyle>
        <color>7dff0000</color>
      </PolyStyle>
    </Style>
    <Placemark>
      <name>Absolute Extruded</name>
      <visibility>1</visibility>
      <styleUrl>#unknown</styleUrl>
      <Polygon>
        <extrude>1</extrude>
        <tessellate>1</tessellate>
        <altitudeMode>absolute</altitudeMode>
        <outerBoundaryIs>
          <LinearRing>
            <coordinates>
              -112.3396586818843,36.14637618647505,1784
              -112.3380597654315,36.14531751871353,1784
              -112.3368254237788,36.14659596244607,1784
              -112.3384555043203,36.14762621763982,1784
              -112.3396586818843,36.14637618647505,1784
            </coordinates>
          </LinearRing>
        </outerBoundaryIs>
      </Polygon>
    </Placemark>
    <Placemark>
      <name>Absolute Extruded</name>
      <visibility>1</visibility>
      <styleUrl>#transBluePoly</styleUrl>
      <Polygon>
        <extrude>1</extrude>
        <tessellate>1</tessellate>
        <altitudeMode>absolute</altitudeMode>
        <altitude>2000</altitude>
        <outerBoundaryIs>
          <LinearRing>
            <coordinates>
              -112.3396586818843,36.14637618647505,2000
              -112.3380597654315,36.14531751871353,2000
              -112.3368254237788,36.14659596244607,2000
              -112.3384555043203,36.14762621763982,2000
              -112.3396586818843,36.14637618647505,2000
            </coordinates>
          </LinearRing>
        </outerBoundaryIs>
      </Polygon>
    </Placemark>
  </Document>
</kml>

Here's a screenshot (using the KML interactive sampler):

I'm dissatisfied with the polygons intersecting color-wise.

(There is also the limitation of extruding up from the surface -- my workaround is to set an altitude above the surface of the earth at which the "extrusions" will start, and then extrude downward from that altitude, so the rock layers will float above the surface.)

If there are go-to ways to do this in ArcGIS or some other package I haven't considered, please disregard my KML approach and suggest the easier way!


texts: Install GDAL extension on GeoServer 2.9-SNAPSHOT

texts:
I have GeoServer on Ubuntu 14 as part of a Geonode 2.6 install.  I am having trouble installing the GDAL extensions.  The version of GeoServer says 2.9-SNAPSHOT.

I have tried to install the GDAL extensions following the instructions here - http://docs.geoserver.org/2.9.3/user/data/raster/gdal.html

I used the GeoServer 2.9.4 GDAL extensions. After restarting GeoServer, the additional coverage formats aren't available.

gdalinfo --version returns the correct version, and
gdalinfo --formats returns the list of supported formats.

I can't find anything obvious in the log file, there is a reference to -
 ArcSDE Java API seems to not be on your classpath. Please verify that all needed jars are. ArcSDE data stores will not be available.

and then

2017-06-27 03:31:35,393 DEBUG [gdal.envihdr] - EnviHdrFormatFactory is not availaible.
2017-06-27 03:31:35,393 DEBUG [gdal.jp2mrisd] - JP2MrSIDFormatFactory is not availaible.
2017-06-27 03:31:35,393 DEBUG [gdal.jp2ecw] - JP2ECWFormatFactory is not availaible.
2017-06-27 03:31:35,393 DEBUG [gdal.mrsid] - MrSIDFormatFactory is not availaible.  ...........etc

I have updated the \etc\environment file to include the path, still nothing.


texts: Function get_partition(geometry, unknown) does not exist - Nominatim tiger data import

texts:
Ive install Nominatim with postgres9.3 and postgis 2.1, imported osm data ok and now trying to import tiger data (tiger2013 edges) for the house numbers.

The import looks to go ok, but there is no data, and the log for postgresql has:

UTC ERROR: function get_partition(geometry, unknown) does not exist at character 8

UTC HINT: No function matches the given name and argument types. You might need to add explicit type casts

UTC QUERY: SELECT get_partition(place_centroid, 'us')

I looked at the source code and in functions.sql it has

CREATE OR REPLACE FUNCTION get_partition(in_country_code VARCHAR(10)) RETURNS INTEGER.

EDIT
_
I asked this over on  https://help.openstreetmap.org/questions/34315/issue-importing-tiger-data-into-nominatim and it turns out there was a bug in nominatim, I should have used the stable version rather than grabbing the latest copy from git


texts:
I want to convert this XML file to a polyline shapefile.

<?xml version="1.0" encoding="utf-8" standalone="yes"?>
<ci nc="BMSKDMA10" nl="BMskDMA10" mt="6" ct="DM" dtt="17342" dtu="15/01/2016 09:47:29">

</ci>

I tried this code with pyshp library, it works as points when I put:

w = shapefile.Writer(shapefile.POINT)

and when I put :

w = shapefile.Writer(shapefile.POLYLINE)

I do not get a polyline

from xml.etree import ElementTree
import shapefile
import os

xml_file = 'BMSKDMA10.xml'
shape_file = 'BMSKDMA10.shp'
projection = 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'

tree = ElementTree.parse(xml_file)
w = shapefile.Writer(shapefile.POLYLINE)

# create fields
w.field('la', fieldType = 'N', size = 200)
w.field('lo', fieldType = 'N', size = 200)
w.field('tr', fieldType = 'C', size = 200)
w.field('nr', fieldType = 'C', size = 200)
w.field('ce', fieldType = 'C', size = 200)
w.field('er', fieldType = 'N', size = 200)
w.field('mt', fieldType = 'N', size = 200)
w.field('cm', fieldType = 'N', size = 200)
w.field('id', fieldType = 'N', size = 200)
w.field('dt', fieldType = 'N', size = 200)

root = tree.getroot()
ci = root.getchildren()

for tr in ci:

    part = []
    part.append([float(tr.get('lo')), float(tr.get('la'))])
    w.poly(parts=[part])

    # copy attributes
    w.record(float(tr.get('la')),float(tr.get('lo')), tr.get('tr'),tr.get('nr'),tr.get('ce'),int(tr.get('er')),int(tr.get('mt')),int(tr.get('cm')),int(tr.get('id')),int(tr.get('dt')))
w.save(shape_file)

# create the PRJ file
with open(os.path.splitext(shape_file)[0] + os.extsep + 'prj', 'w') as prj:prj.write(projection)

Can you help me to resolve this problem?


texts:
I am trying to convert a Google Map KML file to a csv, to map the longitude and latitude points on Leaflet. But the coordinates keep mis-projecting and some end up off the land mass.

The steps I have done are:
- Download KML file
- Import into a GIS package (both MapInfo 12 and QGIS), and used the projections of WGS84 EPSG 3857, WGS84 EPSG 4326 and Google Mercator EPSG 900913, to create long/lat points
- Save this table as csv
- Upload to Leaflet using the Recline.js in DKAN to allow users to preview the map data.

Screenshot 1 is the original google map, and then screenshot 2 is the csv on Leaflet

Have looked around the rest of the forum and cannot find an answer. Any help?


texts: Getting a reference of an already running session of ArcMap with ArcObjects/C#

texts:
I'm new to both C# and ArcObjects and I want to get the layers on a currently opened document, which I've done via:

IApplication app = ???;
IMxDocument doc = (IMxDocument) app.Document;
IEnumLayer layers = doc.FocusMap.Layers;

I have figured out how to get IApplication app if the program opens ArcMap itself[1] or if it's an ArcMap AddIn[2].

I've been trying to figure out how to reference an already running ArcMap document with a stand alone program but have come up empty handed. Is this possible?

[1] https://desktop.arcgis.com/en/arcobjects/latest/net/56eba175-54c1-46b4-801a-fc968e582581.htm
[2] IApplication app = ArcMap.Application;


texts:
I'd like to fix the maximal area, in which the user can fly around. Now its possible for the user, to pan everywhere he wants to go. For e.g. to XMin: 216.2.

My map has got these "borders":

  XMin: 674580.199070543
  YMin: 241290.573962568
  XMax: 691333.99095178
  YMax: 254627.474078552

How can I specify that the user can only pan from XMin to XMax and from YMin to YMax.
I'm using ESRI API for WPF and C# in Visual Studio 2010.


texts: Multiple conditional statements within a Con Tool

texts:
Well I want to apply multiple conditional statements within a Conditional tool. My aim is to loop through a list of rasters and on the basis of a specified range assign each range with a constant 1, 2, 3, 4 or 5 value so that all the final rasters have pixels having a DN of 1,2,3,4, or 5. I copied the template from the following Desktop Help website and adjusted to my requirements.

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Conditional_evaluation_with_Con/009z00000008000000/

My code is as follows:

import arcpy
from arcpy import env
# Set the current workspace
env.workspace = r"C:\Experiments\Test_Jan03-11_A\Masked"
# Get a list of Rasters
rasters = arcpy.ListRasters("*", "ALL")
for raster in rasters:
   output = arcpy.sa.Con(raster <= -10, 1, Con(raster <= 73, 2, Con(raster <= 196, 3,     Con(raster <= 403, 4, Con(raster <= 600, 5, 0)))))
output.save("C:/zE/ResearchWork/UnderWorks/Experiments/Test_Jan03-11_A/Masked/Classify")
Print = "Done"
Runtime error <type 'exceptions.NameError'>: name 'Con' is not defined

The error indicates that I am making a mistake with the conditional statement that I can't identify.

The following code has worked. But with some issues stated below

import arcpy
from arcpy import env
arcpy.env.extent = "MAXOF"
# check out spatial analyst extension
if arcpy.CheckExtension('Spatial') == 'Available':
     arcpy.CheckOutExtension('Spatial')
else:
    arcpy.AddMessage("Error: Couldn't get Spatial Analyst extenstion, exiting")
sys.exit(1)
from arcpy.sa import *
# Set the current workspace
env.workspace = r"C:\zE\ResearchWork\UnderWorks\Experiments\Test_Jan03-11_A\Masked\CMasked"
# set cell size
arcpy.env.cellSize = r"C:\zE\ResearchWork\UnderWorks\Experiments\Test_Jan03-11_A\Masked\CMasked\GW1AM2_20130100_01M_EQMA_L3SGSMCLA1110100"
# Get a list of Rasters
rasters = arcpy.ListRasters("*", "ALL")
for raster in rasters:
    inRas = Raster(raster)
for raster in rasters:
    output = arcpy.sa.Con(inRas <= -10,1,Con(inRas <= 73,2,Con(inRas <= 196,3,Con(inRas <= 403,4,Con(inRas <= 600,5,0)))))
output.save("C:/zE/ResearchWork/UnderWorks/Experiments/Test_Jan03-11_A/Masked/CMasked/MaskedClassify")
Print = "Done"

The problem is that the result is only one raster. I was expecting each of the input rasters to be classified and saved separately. Also it is giving me the following error for saving data.

Runtime error <type 'exceptions.RuntimeError'>: ERROR 010240: Could not save raster dataset to C:/zE/ResearchWork/UnderWorks/Experiments/Test_Jan03-11_A/Masked/CMasked/MaskedClassify with output format GRID.


texts:
I am very confused! I have tried this afternoon to create a polygon of a set size, namely 30m x 1.8m. I need 33 of these set over a grid (which I have created on a vector layer).

I've read everything I can on the forum and tried and failed to use several plugins to do it (e.g. ArcheoCAD, CadDigitize, CADTools) and and now very confused! Some plugins that were mentioned appear to be old and I can't find them (ImprovedPoly? and CADInput). I am running Pisa on a Windows 8 system.

I do not have co-ordinates for the polygons as they are trenches which I simply need to lay on my vector layer a N-S, E-W basis and I was going to use my grid to help me do this.

I have to draw trench plans like this a lot and I am sure there is a very simple way I can do this in QGIS - I have done similar things in CAD.


texts: Creating a 3D model of a pond using GPS coordinates and associated depth measurements

texts: There is a pond which I am trying to create a 3D model of to view the depth. We have measured the depth every few meters along straight line transects all across the pond, and have GPS coordinates for each of these depths. I have all of the data in excel, and have tried using qgis to map the points, but using the .csv data I was only able to create a 2D map of the pond without the depths being entered. How do I add a z-value and create a 3D model of this pond? Coordinates are in degrees, depth is in meters.


texts: ST_NumInteriorRings produces no result in a geometry that has inner rings?

texts:
I am working in Postgres 9.5. This is my table:

 gid    | integer                 | not null default nextval('mytable_gid_seq'::regclass)
 geom   | geometry                |
 county | character varying(1000) |

The geom field was created by running ST_Union over a MultiPolygon field in another table.

If I run SELECT geom FROM mytable WHERE county='Manchester' and then export the results to GeoJSON and open them in QGIS, the resulting polygons look like this.

Now I want to find the number of interior rings in the geometry for Manchester:

SELECT ST_NumInteriorRings(geom) FROM mytable WHERE county='Manchester';

However, it produces an empty result. I don't understand this, since I can clearly see holes in QGIS.

I've also tried ST_NumInteriorRings(ST_Multi(geom)) but that doesn't produce a result either.

What am I doing wrong? (What I would ultimately like to do, in case it helps, is find the area of the internal holes.)


texts: Converting X,Y coordinates to lat/long using pyproj and Proj.4 returns the wrong coordinates

texts:
I'm writing a python script that reads multiple XML files containing x and y coordinates and combines them all into a single csv file.  Latitude and Longitude are required fields in the csv, but I am having difficulty converting the x,y coordinates in Ohio North State Plane usFt to WGS84.

>>> p = Proj(r'+proj=lcc +lat_1=41.7 +lat_2=40.43333333333333 +lat_0=39.66666666666666 +lon_0=-82.5 +x_0=600000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=us-ft +no_defs') #Nad83 State Plane Ohio North US Feet Proj object using parameters
>>> p(739400.91,2339327.3,inverse=True)
(-80.138057868777224, 60.278230707978487)
>>> p1 = Proj(init="epsg:3734") #Nad83 State Plane Ohio North US Feet Proj object using EPSG code
>>> p1(739400.91,2339327.3,inverse=True)
(-80.138057868777224, 60.278230707978487)

Both methods above return the same result, however this lat long is somewhere in Hudson Bay.  When I plot the coordinates in ArcMap, the correct lat long is: -81.142311,41.688205.

*Notice all lat longs are provided long,lat as this is the order Proj uses

Does anyone know why I would be getting the wrong coordinates from Proj.4 and pyproj?


texts:
My Atlas-layer ist "xblattschitt-fuer-rodung-f". In the column "FOTO1" is the name of the jgp. My source is in a network at "z:\01\P156. In the next step I insert a html-Label

<style>
* {margin: 0px; padding: 0px}
</style>
<img src="file:///z/01/P156/[% "FOTO1" %]"
style="width: auto; height: 100%; display: block; margin: 0px 0px 0px auto;" />

Without success


texts: Converting triangle mesh polygon into square cells

texts:
I have a triangle mesh polygon feature class.

I would like to convert it to raster-like cells but still as a polygon feature class using ArcGIS Pro.

This needs to be at 1m resolution

It is important to retain the attributes of the triangles.

I tried doing it like this:

Polygon to raster tool on the triangle mesh
Raster To Point
Create fishnet with input as the Raster
Feature to Polygon on fishnet layer using the layer from  Raster to Point as label features

However, this was still running after 24 hours.

Is there a better method to go about this?


texts: Minimum and maximum coordinates of a polygon in Earth engine

texts:
I have a Fusion table having shape-files of various districts of India and I am using this fusion table in Earth Engine code.

var district = ee.FeatureCollection('ft:1PA2zwArj8EsplrX9eMxJ2H_TICyyx855KPnbJhC1','geometry')
    .filter(ee.Filter.eq('name','Begusarai'));

var district_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
    .filterBounds(district)
    .filterDate('2018-03-01','2018-05-01')
    .filter(ee.Filter.lt('CLOUD_COVER',5))
    .sort('CLOUD_COVER')
    .map(maskClouds)
    .median();

Then I perform some classification on the given image and export the classified image to my drive
'input' variable is having that classified image..

Export.image.toDrive({
  image: input.clip(district),
  description: 'begusarai2018',
  scale: 30,
  region: district
});

I now want to get the maximum and minimum coordinates (Latitude and Longitude) of these districts.
How can I achieve this?


texts: Ingesting Project Assets With Service Account Via Manifest (NodeJs)

texts:
I'm using a service account to connect to Earth Engine API and attempting to ingest some assets from Google cloud Storage using a manifest:
uploadManifest: any = {
    id: `project/martynsproject/assets/${item_id]}`,
    tilesets: [{
        sources: [{
            primaryPath: gsFilePath
        }]
    }],
    pyramidingPolicy: "MEAN",
    properties: {
        ...metadata.properties,
        "system:time_start": aquired.valueOf(),
    }
};

ee.data.newTaskId(1, (taskIds) => {
    ee.data.startIngestion(taskIds[0], uploadManifest, async (result, error) => {
        if (error) {
            console.error(`Ingestion request failed`, error);
        } else {
            var assetId = uploadManifest.id;
            var taskId = result.taskId;
            console.log(`Ingestion Task for file [${gsFilePath}] started. Destination [${assetId}]`);
        }
    });
});

The ingestion task is started successfully, however when I later check the status of the task it has failed with an error:

Asset 'projects/martynsproject/assets' does not exist or doesn't allow this operation.

I'm sure this is because the projects/martynsproject/assets folder doesn't exist, but how do I create this and give appropriate permissions on it for myself and the service account to access it?
*Please note the name of our project is not actually martynsproject, I have just changed it for privacy reasons.


texts: GetFeature WFS request with Spatial Filter throws an error

texts:
I'm using the following request in order to get the closest feature from a layer given some specific coordinates:

http://www.juntadeandalucia.es/medioambiente/mapwms/REDIAM_PPHH_2012/wfs?&VERSION=1.0.0&SERVICE=WFS&REQUEST=GetFeature&TYPENAME=MSPF_POLY&Filter=<Filter><DWithin><PropertyName>msGeometry</PropertyName><gml:Point<gml:coordinates>158380,4130457</gml:coordinates></gml:Point><Distance units='m'>10000</Distance></DWithin></Filter>

But the result is:

mapserv(): Web application error. Traditional BROWSE mode requires a TEMPLATE in the WEB section, but none was provided.

I also have tried with the following example of the documentation a the result is the same:

https://demo.mapserver.org/cgi-bin/wfs?&VERSION=1.0.0&SERVICE=WFS
  &REQUEST=GetFeature&TYPENAME=cities&Filter=<Filter>
  <DWithin><PropertyName>Geometry</PropertyName><gml:Point>
  <gml:coordinates>135.500000,34.666667</gml:coordinates>
  </gml:Point><Distance units='m'>10000</Distance></DWithin></Filter>

I would appreciate some help.


texts: How to create a point along a line given a distance from another point

texts:
I'm trying to figure out a point(D) on B-C that is a specific distance from A.
I really don't have any idea on how to proceed with the ellipsoidal math so I figured I'd ask here. How would I figure out the longitude, latitude of D which is 5mi(26400ft) away from A, on the line B-C? The bearing of A-D would be just as useful as i could figure out the point 5mi away from A with the bearing.

I wouldn't mind if the formula lined up with the Vincenty model of the earth.

Using the Vincenty formula's I've already calculated the following.

Points (longitude, latitude)
A: -86.32899412316736,41.12507719625815
B: -86.29237110757369,41.16845596588609
C: -86.19689279712584,41.16923568742915

A-B Distance: 18750ft
B-C Distance: 26289ft
A-B Initial Bearing: 32.5331889647
B-C Initial Bearing: 89.3493884487

Thanks for the help!


texts: Algorithm to draw a perpendicular line on Google earth (KML)

texts:
I have two points say A (lat,lon ) and B(lat,lon) and I want to draw a perpendicular bisector line to AB such that any point (C) on the perpendicular line satisfies AC=BC.

All I want is a function which takes A (lat,lon) B (lat,lon) and distance from the mid point of A and B and returns lat,lon (C) such AB = AC.

The end goal is to draw various circle with AB as chord.


texts:
I'm trying to draw a combination of layers with Mapserver. One of the layer takes a shapefile as input and the other one is based on GRIB data.

The issue is that the GRIB layer (in green) is drawn from 0 deg to 360 deg. The shapelayers draws from -180 to 180. And indeed, when I examine the GRIB file in Panoply, the meta data of sample shows, among other things, this:

netcdf file:/D:/gribs/sample.grb {
  dimensions:
    lon = 360;
    lat = 181;
    isobaric = 4;
    time = 3;
  variables:
    float lat(lat=181);
      :units = "degrees_north";

    float lon(lon=360);
      :units = "degrees_east";

  ...
  other variables
  ...
  }

Which leads me to conclude that the coordinate system runs from 0 to 360 degrees in the east ward direction and from 0 to 180 in the north direction. So Mapserver is correct when it draws this image:

So my question is how can I wrap the GRIB image so that the right half is placed over the left half of the earth?

The url with parameters i use:

http://mapserver-dev.de.meteogroup.net/cgi-bin/mapserv?map=/var/www/mapserver/routeguard/poc_test.map&mode=map&layer=grib_sample&layer=world_outline_red

My mapfile:

    MAP
    EXTENT      -180 -90 360 90 #max x should be only 180, but then only half of the grib layer is drawn
    CONFIG      "CPL_DEBUG" "ON"
    CONFIG      "MS_ERRORFILE" "/var/www/mapserver/routeguard/poc.log"
    SYMBOLSET   "symbols.sym"
    FONTSET     "../fonts/fonts.txt"
    IMAGETYPE   PNG8
    IMAGECOLOR  255 255 255
    SIZE        1600 900

# ------------------------------------------------------------------------------
# the Map Projection
# ------------------------------------------------------------------------------
    PROJECTION
        #which to use????
        "init=epsg:4326" #54004 #4326
        "units=m"
    END

# ------------------------------------------------------------------------------
# the Web object
# ------------------------------------------------------------------------------

    WEB
        IMAGEURL    '/mapserver/routeguard/' #do not include the the var/www
        TEMPLATE    'poc_template.html'
    END

#--------------------------------------------------------------------------------
# the country mask
#--------------------------------------------------------------------------------
    #red world, defined in shp file
    LAYER
        NAME "world_outline_red"
        TYPE POLYGON
        DATA "/var/www/mapserver/data/routeguard/static/world_outline.shp"
        CLASS
            COLOR 255 0 0
        END

        PROJECTION
            "init=epsg:4326"
        END
    END

#----------------------------------------------------------------------------
#grib example
#----------------------------------------------------------------------------
    LAYER
        NAME "grib_sample"
        TYPE LINE
        CONNECTIONTYPE CONTOUR
        DATA "/var/www/mapserver/data/grib/routeguard/sample.grb"

        PROCESSING "BANDS=2"
        CLASS
            STYLE
                WIDTH 1
                COLOR 0 100 0
            END
        END

        PROJECTION
            "init=epsg:4326"
        END
        METADATA
            #"wms_srs"       "epsg:54004"
        END

    END

END #end map file


texts: Getting the outer boundary of a group of polygons

texts:
I have an issue where I have feature classes of groups of zip codes in different sectors of the US that I have created. I want to be able to define the outer boundary of these groups so I can create a single shape file (one polygon). I was going to use this new polygon to clip the original file.

Is there a way to take many polygons and create a file that only shows the outer extent of the group? And how would I go about such an endeavor?


texts: Conditional formatting attribute table in QGIS print composer

texts: Can you make the same conditional formatting on an attribute table added to the printcomposer as you can do within the attribute table? (with QGIS 2.16.3)


texts: How can GeoServer transform a 'no RHR layer' into a RHR layer?

texts:
We are using GeoServer 2.11.2 and Mapbox.

We have many layers that do not follow the rule of the right hand and we need to serve them to other application. The problem is that we have to show in a map with [Mapbox][1] using the function addLayer and this function only shows layers that follow RHR.

Is there any possibility of serving an RHR directly with GeoServer from a no RHR layer? Or is there other function in Mapbox that allow showing these kinds of layers?

I'm using GeoJSON as output format


texts: Calculate the area of a map using gdalinfo

texts:
I'm using gdalinfo to get some information out of a map, and want to use it to calculate the area represented by the map. I'm currently doing it (when the units are given in metres), by multiplying size[0]*size[1]*geoTransform[1]*-geoTransform[5]/1,000,000 to get the area in square kilometres.

I've manually confirmed many of these areas by entering the wgs84extent into geojson. However, this usually only works when the skew and rotation values (geoTransform[2] and geoTransform[4]) are 0. If they are non-zero, my area calculation is wrong. I was wondering what equation I should be using to consider these values and get an accurate area?

This is an example. The units are again in metres here. My calculations get me 0.025 square kilometres, while geojson gives me 16.3.

  "size":[
    2565,
    2052
  ]

  "coordinateSystem":{
    "wkt":"PROJCS[\"NAD_1983_UTM_Zone_11N\",\n    GEOGCS[\"GCS_North_American_1983\",\n        DATUM[\"North_American_Datum_1983\",\n            SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],\n        PRIMEM[\"Greenwich\",0.0],\n        UNIT[\"Degree\",0.0174532925199433]],\n    PROJECTION[\"Transverse_Mercator\"],\n    PARAMETER[\"False_Easting\",500000.0],\n    PARAMETER[\"False_Northing\",0.0],\n    PARAMETER[\"Central_Meridian\",-117.0],\n    PARAMETER[\"Scale_Factor\",0.9996],\n    PARAMETER[\"Latitude_Of_Origin\",0.0],\n    UNIT[\"Meter\",1.0]]"
  }

  "geoTransform":[
    569166.1990436509950086,
    0.070050043573309,
    -1.7578549919737869,
    4986248.3334737550467253,
    -1.7582010122098635,
    -0.0705420294500134
  ]

  "wgs84Extent":{
    "type":"Polygon",
    "coordinates":[
      [
        [
          -116.1220248,
          45.0263128
        ],
        [
          -116.1678264,
          45.0253527
        ],
        [
          -116.166135,
          44.9847436
        ],
        [
          -116.1203658,
          44.9857033
        ],
        [
          -116.1220248,
          45.0263128
        ]
      ]
    ]
  }


texts: gdalwarp -r max and -r average giving the same (incorrect) result

texts:
I am resampling a raster from ~500m to ~1km in using two different methods: max and average. However, both results are the exact same and not what I'd expect for either.

gdalwarp -ot Byte -tr 926.625432000000046 926.625432000000046 -r max input.tif output-max.tif

gdalwarp -ot Byte -tr 926.625432000000046 926.625432000000046 -r average input.tif output-ave.tif

input.tif is a raster consisting of 0's, 100's, and NoData. Its resolution is 463.312716000000023, exactly half of the output resolution. Therefore I would expect the results of the max resampling to be 0 or 100, and the results of the average resampling to be 0, 25, 50, 75, or 100. However both results are the exact same with values of 0, 25, 33, 50, 67, 75, 100. This is not at all what I would expect. To the naked eye at least the pixels do match up (i.e. four input pixels fit perfectly inside one output pixel) so I'm not seeing any reason for this behavior. Any ideas?


texts:
I have a massive Dataset from OpenStreetMap and imported it into a PostgreSQL Database using ogr2ogr. I got several warnings during this process because of multipolygons with too little vertices. Now I can’t transform the data in the Database using PostGIS, because of this invalid features. It is no Problem to losse the invalid features but I don't know how to delete them or correct them. I tried:
ALTER TABLE multipolygons ALTER COLUMN wkb_geometry TYPE geometry USING ST_Transform(ST_MakeValid(ST_SetSRID(wkb_geometry,4326)),3857);
But I got the error Message:
Polygon must have at least four points in each ring
Any recommendations?
I tried:
Alter table multipolygons ADD COLUMN wkb_geom_valid geometry; Update multipolygons SET wkb_geom_valid = ST_MakeValid(wkb_geometry);
delete from multipolygons where ST_NPoints(wkb_geometry)<4 ;
select * from multipolygons where ST_isvalid(wkb_geometry)=TRUE ;
All without success.


texts: Data defined label position for more than one scale in one layer

texts:
I'm looking for a better labeling-solution for my projects. It seems to be difficult - but maybe I'm just blind?

My projects:
I build  printed bike-hike-books. The main scale is 1:75k and there are a city scale of 10k and a overview scale of 200k. Because I like to have a "cartographic" quality of labels placements, I touch every single label and move it to it's best position. I have to do this for every single scale so my labeling columns for 3 scales look like this ('scale 1' stands for 10k 'scale 2' for 75k  and 'scale 3' for 200k):

But before I try to explain in detail, I build this simple Example-Project where you can see and test my (old) labeling solution for 3 scales: http://www.map-site.de/tmp/LabelingExample.zip

Well, it is working but:

You need labeling columns for every scale
You need a layer for every scale

I'm looking for a much better solution without duplicate layers and tons of columns. I count on your creative minds!

EDIT: I added this two very dense map as examples:

There is only one single possible position for "St. Paulus Dom" and it's not any quadrant:

Every label has to be moved manually. Some of the labels are not in a quadrant.


texts: Linux (unix) bash, import multiple fgdb to PostGIS using find & ogr2ogr command

texts:
I have multiple file geodatabases in different folders and this is the command I made so far:

find -path "*/V*/*.gdb" -o -path "*/W*/*.gdb" -o -path "*/R*/*.gdb" -o -path "*/WM*/*.gdb" | while read filename;
do ogr2ogr -update -append -progress --config PG_USE_COPY YES --config FGDB_BULK_ROAD YES -f "PostgreSQL" \
PG:"host=localhost port=5432 dbname=postgres user=postgres password=postgres" ${filename} done;```

Here, when I just typed:

find -path "*/V*/*.gdb" -o -path "*/W*/*.gdb" -o -path "*/R*/*.gdb" -o -path "*/WM*/*.gdb" | while read filename; do echo $filename; done;

then I get ./Vblahbah/file.gdb

Since *.gdb is being recognized as folder in Linux, I think using find -path is the best option that I can think of

but when I executed the command I get this error

  Unable to open datasource './RRLM171G0WRR000FZZBB/file.gdb' with the following drivers.

but it does look like the file geodatabase is being imported to PostGIS.

What am I doing wrong?


texts:
I'm working in both ArcMap and Grass

I have a problem with the boundaries or the edges of my vectors. The vectors are a result of object segmentation (called classification) in Grass Gis. In Grass, an images is segmented and then I converted the raster into vector resulting in this:

(this has been edited) And this is great however, the vector layer is framed. The 'frame' isn't a separate thing; those lines are polygon boundaries like any other on the interior.
There is a square around the vector layer but because I have a little over thousand object segmentation rasters (which I'll turn into vectors) I need continuation, the frame must go.

Any thoughts on how to eliminate the frame? I'm open to Grass and ArcMap solutions. Maybe the solution is somewhere earlier in the process (the process being object classification to vector)

Edit:
Here is an example of the intersection of 4 vector layers.


texts: Geographic Center of Earth centric Projections

texts:
This is more a question born out of curiosity than actual needs, but this seems like the right place for it.

Based on What is Centroid of all lands of Earth? I was wondering what world maps centered on the GC would look like. Has somebody played around with it and can link their results here(or even place it on wikipedia in this article)? I imagine a polar projection or a cylindrical projection with the GC as coordinate origin.


texts: Function for hyperspectral rasters with 3 variables unexpected error

texts:
I am writing a script to extract vegetation indices from hyperspectral datacubes. One of the indices is the Triangular Vegetation Index for which I wrote the following function:

    TVI2 <- function(h,m,l)
{
  (1.5*(2.5*(h-m)-1.3*(h-l))/sqrt((2*h+1)^2-(6*h-5*sqrt(m)-0.5)))
}

The variables stand for different wavelengths in the forumula. I have a rasterbrick object containing all 8 rasters and 101 bands per raster. To call a specific raster and wavelength I use:

raster(r[[1]], layer=1)

And this works just fine for calculating an index with just 2 variables, like the NDVI function:

ndvi1 <- NDVI(raster(r[[1]], layer=41),(raster(r[[1]], layer=71)))

However if I want to use this same set up for the TVI2 function:

tvi1 <- TVI2(raster(r[[1]], layer=71),(raster(r[[1]], layer=45),(raster(r[[1]], layer=21))))

I keep getting the following error:

Error: unexpected ',' in "tvi1 <- TVI2(raster(r[[1]], layer=71),(raster(r[[1]], layer=45),"

I cannot seem to figure out why. If I don' t use the raster call as the last variable, but just a number, the function runs fine.

Any ideas on how to solve this?


texts: Static image layer combined with vector layer in OpenLayers 6

texts:
I'm trying to have an image layer combined with a vector layer but when I add projection to View the vectors disappear. Am I doing something wrong?

Here is an example of what I'm trying to achieve: https://stackblitz.com/edit/react-ol6

I want these two layers below to be together in View.

Vector Layer

Image Layer

import React, { Component } from "react";
import ReactDOM from "react-dom"

import "./style.css";
import "ol/ol.css";
import { Map, View } from "ol";
import { getCenter } from "ol/extent";
import ImageLayer from "ol/layer/Image";

import Feature from "ol/Feature";
import Point from "ol/geom/Point";
import VectorLayer from "ol/layer/Vector";
import TileLayer from "ol/layer/Tile";
import VectorSource from "ol/source/Vector";
import Projection from "ol/proj/Projection";
import Static from "ol/source/ImageStatic";

import { Fill, RegularShape, Stroke, Style } from "ol/style";

var stroke = new Stroke({ color: "black", width: 2 });

var styles = {
    square: new Style({
        image: new RegularShape({
            stroke: stroke,
            points: 4,
            radius: 10,
            angle: Math.PI / 4,
        }),
    }),
};

var extent = [0, 0, 2048, 2048];

var styleKeys = "square";
var count = 50;
var features = new Array(count);
var e = 450000;
for (var i = 0; i < count; ++i) {
    var coordinates = [2 * e * Math.random() - e, 2 * e * Math.random() - e];
    features[i] = new Feature(new Point(coordinates));
    features[i].setStyle(styles[styleKeys]);
}

var source = new VectorSource({});

var vectorLayer = new VectorLayer({
    source: source,
});

var imageLayer = new ImageLayer({
    source: new Static({
        url: "https://imgs.xkcd.com/comics/online_communities.png",
        projection: new Projection({
            units: "pixels",
            extent: extent,
        }),
        imageExtent: extent,
    }),
});
vectorLayer.setZIndex(1);
imageLayer.setZIndex(0);

var olMap = new Map({
    target: null,
    layers: [imageLayer, vectorLayer],

    view: new View({
    //if i add projection the vector layer doesnt render

     projection: new Projection({
           units: "pixels",
              extent: extent,
           }),
        center: getCenter(extent),
    center: [0,0],
        zoom: 2,
    }),
});

class App extends Component {
    componentDidMount() {
    olMap.setTarget("olmap");
        vectorLayer.getSource().addFeatures(features);
    }

    render() {
        return (

        );
    }
}


texts:
I have got a .ecw-file that I want to somehow import into QGIS 1.8.0.

Is there a plugin that I can install that can help me do this or is there another way?

The .ecw-file contains raster data.


texts: Transforming an image (TIFF) into an array

texts:
I found the following code for crop classification and I am trying to understand it.

# We have 10 images and 5 bands for each image so we’ll have a total of 50 columns in our data frame
# we’ll sample the 100,000 pixels for each class to avoid overfitting

temp = pd.DataFrame()
final = pd.DataFrame()
for c in top_classes:

    # print('\n Reading for class: {} \n'.format(c))

    #Read image ``img`` and return ``np.array`` of image values Image will be (nband, nrow, ncol)

     for img in images:

        train_ds = gdal.Open(img, gdal.GA_ReadOnly)

        nrow, ncol, nband = train_ds.RasterXSize, train_ds.RasterYSize, train_ds.RasterCount
        dtype = gdal_array.GDALTypeCodeToNumericTypeCode(train_ds.GetRasterBand(1).DataType)

        #print('Image {}:  width= {}, height= {}, number of bands= {}'.format(img, nrow, ncol, nband))

        img_b1 = np.zeros((ncol,nrow,nband),dtype=dtype)

        for b in range(nband):
            img_b1[:, :, b] = train_ds.GetRasterBand(b + 1).ReadAsArray()
        #print(img_b1)

        Xt = img_b1[roi==c, :]

        Xt1 = pd.DataFrame(Xt)
        Xt1.head()

       # Xt2 = Xt1.sample(n=100000)

       # Xt2.reset_index(drop=True,inplace=True)

       # temp = pd.concat([Xt2,temp],axis=1)

        #temp["class"] = c
        #temp.reset_index(drop=True,inplace=True)

    #final = pd.concat([temp,final],axis=0)
    #final.reset_index(drop=True,inplace=True)

   # gc.collect()

Can someone explain why they looped over the number of bands?

 for b in range(nband):
       img_b1[:, :, b] = train_ds.GetRasterBand(b + 1).ReadAsArray()
 print(img_b1)


texts: ModelBuilder model only works in edit window?

texts:
I have a model that iterates through a single field (such as Name) and selects out each name in that field and ultimately creates a KML layer as the last output of each name in the field by using an in line varible (For Example %Value%). The goal is to have a tool that one can input any shapefile, choose their desired field of choice, and have the tool iterate and create KMLs of each individual value in that chosen field. The issue is right now the tool works fine when I open it in the Edit window and add the data there in the ModelBuilder window but I was hoping to be able to run it like a traditional tool in ArcMap by clicking open and adding the data there. I have set every output from each tool in the model as a parameter so I can control where the temporary files are being stored at every step.

How can I get a successful model that runs in the edit window but not from the open window?


texts: gdalwarp clipping produces a black image

texts:
I have a large radar image (1.1GB) from Sentinel-1 and I need to make it smaller so that I can process it without my memory crashing. Using QGIS I generated the command to mask an image based on a shapefile using Raster > Extraction > Clipper.

The command was gdalwarp -q -cutline my_shapefile -tr 0.000107418644212 0.000107418644212 -of GTiff my_sourceTIF my_destinationTIF.

When I run this command in my OSX terminal it generates a new TIF image. But when I open it in QGIS, the image covers the same area as before except the whole thing is now black.


texts:
I am developing a GIS application like this:

By now my steps are:

I have got a jpg file, 5040 x 8179 dimension. In this map, the roads, rivers, towns etc. are marked.
I am using .Net to develop and this requires a shp file to be used properly.
I tried using qGIS to vectorize this jpg and georeferencing it but the shapefile that comes out is like full green color etc..

How to convert my original jpg map data to shp format?

I have very little experience in GIS.


texts: Nearest distance between point and line layers in QGIS

texts:
I have two layers, a line and a point layer.

I need to make a line layer which connects all the points to the nearest line feature. How can I do that?

Is there a plugin available for QGIS?
This is a very important tool which is missing in QGIS.

ArcView has this tool: Nearest Features


texts: Can Extent Indicator symbology by automated with ArcPy?

texts:
I have a map series of 30+ data driven pages that I am altering the layout of, trying to do as much of it with ArcPy as possible to save time. I need to change the symbology of the Extent Indicator.

Can this be coded?

Original:

Desired:


texts: Improving speed while adding QgsVectorLayer

texts:
In my QGIS plugin I'm adding some different layers to the QGIS canvas with a uri connection to my Postgres database, adding one layer takes no time, hence adding 20 layers takes some time. So I tried to understand what consumed the time by adding some t=time.time(), this is part of my code:

host = self.dbhost
port = self.dbport
dbname = self.dbname
username = self.dbuser
password = self.dbpass
t0 = time.time()
## Adds a PostGIS table to the map
uri = QgsDataSourceUri()
uri.setConnection(str(host), str(port), str(dbname), str(username), str(password))
uri.setDataSource(str(schema), str(table), str(geom_col), filter_text, 'field_row_id')
uri.setKeyColumn('field_row_id')
uri.setSrid('4326')
t1 = time.time()
vlayer = QgsVectorLayer(uri.uri(), str(extra_name) + str(table), 'postgres')
t2 = time.time()
if not vlayer.isValid():
    QMessageBox.information(None, 'Error', 'Layer not loaded correctly, Connection:\n' + str(uri.uri()))
print('DBfunc: setting uri con {u} s, creating layer from uri {s}, checking it {t1}'.format(u=round(t1-t0,3),s=round(t2-t1,3),t1=round(time.time()-t2,2)))

So the log clearly showed that setting the Uri connection (the code until t1 took no time (less than 0.001 s), checking if the layer was valid took no time (less than 0.001 s), hence the code line creating the QgsVectorLayer took approximate 0.5 sec to execute. Is here any improvements to be made? (using QGIS3)


texts: Add multiple strokes to Leaflet geojson feature path

texts:
Wondering if it's possible to add multiple strokes to a polygon or line feature (the same thing could apply to multiple fills for polygons too).

The leaflet documentation on the L.Path style options doest seem to suggest there is a way.  For example, styling the stroke once is straightforward:

var geoJsonFeature = L.geoJson(data, {
  style: {
    color: '#ff9700',
    weight: 4
  }
}).addTo(map);

... but it doesn't appear to allow adding a second, thin stroke on top (e.g. with a thinner weight, so the two are strokes are layered one on top of the other), w/o, I guess, loading the geojson twice.

Is it possible to hack this somehow?  Or might it be possible with another library like d3.js?


texts:
Is it possible to write sql files containing all the CREATE TABLE... and INSERT... statements etc. with FME instead of directly execute them against the database (here: PostGIS)?

Background is that I have a PostGIS DB on a remote server that I am not able to access with FME Workbench directly, but I can transfer i.e. textfiles (like *.sql) to it in order to build the databases.


texts: ArcGIS online widgets - filter by attribute

texts: I started to use ArcGIS online. I created my map and I published it as a web app. I would like to ask if there is any possibility to allow users to filter by attributes using widgets or so on? I can't find this option.


texts:
How do I set mouse wheel zooms in ArcMap 10.1?
We want zooming in ArcMap to work like it does in OpenLayers.


texts:
I'm trying to add edit(namely Edit.MOVE, Edit.SCALE, Edit.ROTATE, Edit.EDIT_VERTICES, and Edit.EDIT_TEXT) to the draw widget inside CMV. The steps I have taken thus far are:

1- I've defined the edit namespace:
'esri/toolbars/edit',

2- I've added "Edit" inside the function declaration.

3- Created an editTools variable:
drawToolbar: null,
editTools: null,

4- Added the following code inside the post create method to respond to different types of graphics:
this.editTools = new Edit(this.map);

        this.pointGraphics.on('click', function(evt)
        {
            console.log("CLICKED ON POINT GRAPHIC");

            event.stop(evt);
            this.activateEditTools(evt.graphic);
        });

        this.polylineGraphics.on('click', function(evt)
        {
            console.log("CLICKED ON POLYLINE GRAPHIC");

            event.stop(evt);
            this.activateEditTools(evt.graphic);
        });

        this.polygonGraphics.on('click', function(evt)
        {
            console.log("CLICKED ON POLYGON GRAPHIC");

            event.stop(evt);
            this.activateEditTools(evt.graphic);
        });

        this.map.on('click', function(evt)
        {
            console.log("CLICKED ON MAP");

            this.editTools.deactivate();
        });

5- The activate edit tools method I'm calling above:
        activateEditTools: function(graphic)
        {
            alert("inside activate edit with graphic: " + graphic);
            this.editTools.activate(Edit.MOVE | Edit.SCALE | Edit.ROTATE | Edit.EDIT_VERTICES | Edit.EDIT_TEXT, graphic);
        },

Yet to no avail, any input would be helpful.


texts: Fixing Model with Iterate selected feature?

texts:
I have a Shapefile with point features representing pour points. I need to delineate their watersheds individually, because many of them overlap.

I'm using the following model to select a point from the Shapefile, export it in a new layer, delineate it's watershed with the Watershed tool, and then repeat for the next point. I want each new layer (of point and watershed) to have the FID as the file name, so I'm using %FID% as the output.

It keeps failing. What's wrong with my model? I'm relatively new to ModelBuilder.


texts: Calculate Volume of Thickness/Difference Raster

texts:
I have two LiDAR datasets, and I would like to know the volume of their difference. I used raster calculator to produce a "difference" raster. What I need now is a way to be able to calculate the total volume of that raster. I don't think Surface Volume, nor Cut Fill, is what I am after as they both use a reference plane or require before/after rasters.

Is there a tool that is well suited to this type of calculation?


texts: Data points Imported into QGIS are in the middle of ocean

texts:
I'm new to QGIS and tried importing some GTFS data in the CSV format (stops.txt). The project CRS, base map CRS and GTFS CRS have all been set to Google Mercator, but the points that are supposed to appear in Boston is now plotted in the middle of the ocean.

Did I miss out on something that messed up the import/plotting?


texts: Finding TOC Layers Type (layer,grouplayer,..) and layerindex of a layer

texts:
How can i find the type of the layers in TOC ?
these types are very important to me :

1.layer (that has a valid datasource)

2.layer (that lose its datasource because of changing the address of its datasource )

3.GroupLayer

And i want to know about the index of the layer that we have with ILayer. A function like this :
Public Function GetLayerIndex(mLayer as ILayer) as long


texts: Seeking example or docs for geotk-client-csw package

texts:
I need a simple but fast CSW client. CatalogServicesClient is created terribly slow. I was about to write my own pool but suddenly found org.geotoolkit.csw.CSWClientFactory not mentioned in anywhere. It has open() with sematics that makes no sense to me.

Any directions please?


texts: v.dissolve renames the attribute table column and does not produce a dbf

texts:
When I run v.dissolve on my vector...

v.dissolve --o input=pass2@mapset1 layer=1 column=fttcp output=pass3

The dissolving appears to work but the column name in the attribute table has changed from 'fttcp' to 'cat'.

Also, although the input has a dbf file, the output does not (so I assume it is storing the attribute table in a different format or something?).

Is there a way to tell v.dissolve to retain the original column name or is there a way to rename it?

Also, is there a reason v.dissolve doesnt write a dbf. Is there a way to get it to write the attribute table to a dbf or is there a way to convert the vector so it uses a dbf?

I am using GRASS 6.4.2 on windows 7


texts: I tried to dock the snapping panel by double-clicking on top of the panel but it's still not docking with the main QGIS window.


texts: Export ERA5 Image Collection data from Google Earth Engine

texts:
I'm trying to export the ERA5daily Image Collection data from Google Earth Engine using the GEE API in python. I need the Max Temp, Min Temp and Average Temp for the time between 1979 and 2016 (in the example here just 1 day for convenience), within a specific bbox.

Basically I need the temperature time series (max, min, avg) for each grid within my geometry object, because I want to calculate different indices of temperature variability and extremes afterwards.

Unfortunately I haven't had much success so far. I managed to export single .tif images with the GEE Code Editor in JavaScript but it seems quite unreliable and I would have to press 'Run' thousands of times. I've read that for batch exports the Python API is recommended.

Below is the code I have tried so far, adopted  from GEE Tools

This code works for exporting the images in .tif, however I'm running into several problems. First I don't manage to select the bands. In the GEE code editor this works with

var col = ee.ImageCollection('ECMWF/ERA5/DAILY').select()

But in Python it just doesn't export anything to drive when I use it.

Secondly I can't figure out what the scale parameter does. The default setting is 30, which leads to massive files (and the need to increase maxPixels). When I set it to a higher value however the data stored in the .tif changes. So not really sure what to do about that.

I know there is likely a much better way to extract such large time series to reduce file size. So any help regarding the code adapted from geetools or a better approach for my problem is appreciated.

from geetools import batch
import ee
ee.Initialize()

region = ee.Geometry.Rectangle([79, 25, 89, 31])

col = ee.ImageCollection("ECMWF/ERA5/DAILY")\
    .filterDate('1979-01-02', '1979-01-03')\
    .filterBounds(region)

batch.Export.imagecollection.toDrive(col, "ERA5", region=region)


texts: Merge two rasters with different resolution

texts:
I'm a student and pretty much new with Qgis.

I'm trying to merge two elevation rasters. They have different resolution (one is a 5mx5m DTM and the other is a 1mx1m DTM) and different extension.

I've clipped the two rasters out of two bigger rasters using Clip with mask.

I've tried to resample the lower resolution raster using gal.warp(reproject) but the all system crashes down every time.

Is there a solution or maybe another way to merge them?

Here are the parameters I run the algorithm with:

I also tried to merge them with no resampling using gdal_merge with senseless results. Here's the gdal_merge parameters:

Heading ##


texts:
Is there a way to access the GeoProcessing tools in QGIS, in code? I want to find the Difference between two layers and make a new shapefile with the result.It can be manually done by selecting vector->GeoProcessing Tools->Difference. But how to do it in my code? (I use python)
Please help me to find a solution to this problem


texts:
I've been working on opening up our spatial data from a proprietry system to one that can be read by more products (MapGuide, ESRI, gvSIG, Map3D etc....)

I've decided on using SQL Spatial as we have a server in house available to us and we can link to other systems easily enough through it.

Before I go all out in actually using it as our central datastore is there anything I should be aware of that could greatly decrease performance?


texts: Finding a center of polygons, and its coordinates

texts: So I've got multiple polygons. I need to find the center of that area (not the center of each polygon) and get the coordinates of that point and also get the name of the polygon in which that point is in. I know that with st_centroid I can calculate the center of each polygon, but how do I get a single center point,name of the polygon that it is in and its coordinates?


texts:
I am using OpenLayers and until now I have been using OSM as background and I got points from my database (mySQL) and added them to vector layer, then I could modify them and save it back to the server.

Now I am not using OSM any more, I got data as a shapefile, and this is the background of the map, I am using GeoServer to show the files.

I want to make an option to modify these layers (WMS). I can't insert everything to the database because it's too much and the map will work slowly if I will add the data from the database.

So I need a way to modify the vectors that I am getting from GeoServer and to save the changes, or if its not possible, maybe to do it with KML layer (but this is a bad option).

Is there a way to do it?


texts: ogr2ogr convert DXF to GeoJSON with defined geometry

texts:
I trie to convert a DXF file into a GeoJSON file with:

ogr2ogr -f 'GeoJSON' -sql "SELECT * FROM entities WHERE layer = 'geom_test'" test.json dfxtest.dxf

Is it possible to convert only geometries with for example the type point?

Are there documentations or explanations on how to use ogr2ogr with dxf?


texts: How to work with Baccini et al.'s biomass map from ArcGIS?

texts:
Baccini et al. (2012) has produced one of the best biomass maps available for the world. One can find the map from one of ArcGIS repositories here.

I downloaded the .pkinfo file for that map but I could not work with it. The map is 'unselectable' in ArcMap. The only thing I managed to do was to export it as a .lyr file, but that is not helpful at all.

I am not really a GIS expert so I cannot tell if I missing something really basic about that map but I wonder if there is an easy way to be able to work with that map in ArcMap or even in R? For instance, to be able to crop it and calculate mean values, etc.


texts: How can I activate the Text annotation edit window in QGIS?

texts:
I have a python plugin for QGIS where I manage to show a list annotations auto created by the plugin.

Meanwhile, I can display the annotations in a combobox, on user selection the annotation object is selected on the canvas (and this is done already)

I need to open the annotation text edit window for the selected annotation .. how can I do that .. can anyone guide me to a proper solution ?


texts: Place points to minimize population weighted euclidean distance to census centroids?

texts: In ArcGIS 10.22:  I need to take about 3000 contiguous census blocks and determine where to place 15 points such that the points minimize the population-weighted Euclidean distance to the census block centroids.   I'm trying to think of several ways to do it,  such as creating a surface using population as the z-value.  What would be the most efficient way to calculate this?  Which tools would you use?  (think,  I'm going to drop x piles of corn in the area of varying density such that it would minimize the flying distance for the most crows to fly to those corn piles)


texts:
I'm wondering if it is possible to get zones based on a raster file as you can see in my attached image (hand-made).

Is there a tool or an idea to get a result like that?


texts: Calculating optimal routes while touching all available roads in ArcGIS network?

texts:
For preperation of the winter season we want to calculate the most optimal routes for sprinkling salt on the roads. The analysis knows the following criteria:

vehicles start and stop at a single loading point
all available roads need to be sprinkled with salt
one route can take no longer than a certin time (lets assume 2 hours)
because of the limited load of salt per vehicle, the distance of a route is limited to the amount of salt available. (lets assume 10 km)

Network analyst of ArcGIS (10.0) assumes you have a start- and an endpoint to calculate the route. However, in this case it is not about calculating the fastest route from origin to destination, but about the most optimal routes to cover as much road distance as possible within a limited time frame.

Now we are thinking about calculating midpoints for every road section and use those as destinations to calculate the route.


texts: Geographically closest 2 airports between 2 cities

texts:
I have a database of all airports and need a query that finds the shortest distance between airports in 2 locations (cities for example). Here is my query.

SELECT
  a.city,
  a.iata_faa,
  md.min_distance
FROM
   airports as a, (
        SELECT
          aa.point AS pointa,
          ab.point AS pointb,
          min(st_distance(aa.point,ab.point)) AS min_distance
        FROM
          airports AS aa,
          airports AS ab
        WHERE
          aa.city = 'Miami' AND ab.city = 'Atlanta'
        GROUP BY 1,2

) AS md

WHERE
  st_distance(a.pointa, md.pointb) = md.min_distance
;

I'm closer and tried the <-> operator in the ORDER BY clause but what I'm getting is and ambiguous column reference:

ERROR:  column reference "point" is ambiguous
LINE 23:   st_distance(a.point, md.point) = md.min_distance

Revised query:

SELECT
  a.city,
  b.city,
  md.min_distance
FROM
  airports AS a, airports AS b,
  (
      SELECT
        aa.point,
        ab.point,
        min(st_distance(aa.point,ab.point)) AS min_distance
      FROM
        airports AS aa,
        airports AS ab
      WHERE
        aa.city = 'Miami' AND ab.city = 'Atlanta'
      GROUP BY 1,2

  ) AS md

WHERE
  a.city = 'Miami' AND b.city = 'Atlanta' AND
  st_distance(a.point, md.point) = md.min_distance

GROUP BY 1,2
ORDER BY a.point <-> b.point::geometry limit 1;
;


texts: Generating chainage nodes and labels along line each 10 meters in QGIS

texts:
I have a line defined by a start and end node. I’d like to insert a node every 10m spacing and label it with the chainage (distance along the line).

I thought I could insert nodes every 10m then label it - but I can only find functions to do this for a single node, not multiples.

Can this be completed with the existing functions, or can I leverage some spatialite functions?


texts: Changing Field that Attribute Window orders Layers in using ArcMap

texts:
How do I change the field that orders the attributes?
As you can see, STRUCT_TYP has the same value for each attribute.  I would like to change it to STRUCT_ID.
Is this possible using ArcMap?
I've tried nothing and I'm all out of ideas.


texts: Geopandas and QGIS Calulating Different Polygon Area Values?

texts:
I am just trying to understand why Geopandas and QGIS might calculate areas differently despite using exactly the same shapefiles and exactly the same projections.

Both shapefiles are using EPSG:53032 - Sphere Azimuthal Equidistant. Proj4 string is:

  +proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs

Despite this both have returned similar, but substantially different results.

QGIS Results:

Geopandas Results:

FYI:

QGIS results were calculated using the $area function in the field calculator
Shapefiles are accessible here if you want to play with them.
Geopandas code is below:

import pandas as pd
import geopandas as gpd

# Import Shapefiles from Google Drive
gdf = gpd.read_file('https://drive.google.com/uc?export=download&id=1RnniXEYTJAfazYKQBMIRphj_gloK9uSm')

# Update CRS
newcrs = "+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs"

gdf = gdf.to_crs(newcrs)
gdf["area"] = (((gdf['geometry'].area)/1000)/1000) # Divited by 1000 to convert into metres

gdf.sort_values(['area'], ascending = False)


texts: How to merge 2 polygons from a single layer in QGIS?

texts: I have created buffers around several points, but they come out in patches. Therefore I want to edit and merge them together. How do I do that in QGIS? The red circles shows the areas I want to join.


texts:
I am using offline OSM data for the area of the Czech Republic in my project (freely available for download here). I am using the OSM data in Spatialite + EU-DEM derived contour lines - only as background map.

It is relatively slow even on powerful PCs. I searched for a way to improve the performance and Martin Dobias advised me to optimize the vector data and keep it as vectors.

I plan to reproject the data to EPSG:3857 - to match project CRS and not to complicate switching between offline maps and online via Openlayers.

I would also like to clean the OSM data - remove all unused fields and make it possible to use categorized renderer instead of rule-based renderer.

Is there any user-friendly way to work with the fields similarly to attribute table tools? I found the DB structure quite confusing for me.

Or is it easier to edit directly the OSM format and modify the data structure before conversion to QGIS DB?

Until now I only used CLI commands found in some tutorials to convert OSM to DB.


texts: Intersecting lines and polygons and calculating line length in R?

texts:
I have shapefiles of polygons and lines at my disposal. I would like to intersect polygons and lines with R and thus calculate the length of the lines present in each polygon. You can find the data here "https://www.dropbox.com/s/b6x9gk8bzhan7zm/test.rar?dl=0" .I know how to do it on Arcgis but I would like to have a code on R to automate this because I have to do the same procedure for a lot of data. The code I put in place on R :

setwd("D:/folder/test")

library(sp)
library(maptools)
library(rgeos)
library(rgdal)

Poly <- readOGR("grille.shp",layer="grille")
plot(Poly)

line <- readOGR("line.shp",layer="line")
plot(line)

Inter <- gIntersection(line,Poly,byid=TRUE)

allows me to have an intersection but what I would like to know is how to calculate then the length of the lines present in each polygon ?


texts: Connecting points to form lines, by only connecting points nearest to each other

texts:
I have a layer of points, that I hope to connect to form a polyline(s).  The layer looks like this:

I am hoping the end results traces the network that was originally used to construct the points:

Background:
I am using ArcMap 10.3.  My workflow has been to turn the original polylines into points, after setting a starting point, an interval, and a total stop distance (along each branch in the network).  The I do some processing (data addition) on the points (which is originally a much larger set than shown) outside of arcmap, bring the points back in, create a layer from the initial larger dataset where a certain criteria is set, and then hopefully reverse the points into lines.

If I use the Points to Line tool (Data Management), I get this:

Which is expected since the numbering of the points along so many branches make's it illogical to follow the numbering when creating a polyline.  Is there a way to generate polylines from points when only the nearest points are connected?


texts: QGIS Show several images for a feature (Image Slider widget)

texts:
I have a layer where each feature can have more than one image. Some have 1 photo, some 5 photos, some 15.

I would like to be able to show all available photos of a item in QGIS. The available Photo Widget is only capable os showing one photo.

My photos have the following nomenclature:

00001_01.jpg: first image of feature id 1
00001_02.jpg: second image of feature id 1
00001_03.jpg: third image of feature id 1
00002_01.jpg: first image of feature id 2
00003_01.jpg: first image of feature id 3
00003_02.jpg: second image of feature id 3
.....

I know that there is not a slider widget in QGIS but is there a way to view all the photos for each feature?


texts: Currently I have access to truck fleet data through Telogis Fleet. Telogis provides an XML feed and I would like to know how I can feed that into ArcGIS that is continuously updated?


texts:
We are using a product that makes use of ArcGIS Runtime Engine 10.2. In attempting to find a particular address using the World Geocode Service (ArcGIS Online), I'm receiving nonsensical results:

Looking up: 100 Woodside Cr NW, Airdrie /
Result:
(1) 100 Woodside Cres NW, Airdrie, Alberta, T4B -- 87%
(2) 100 Woodside Cir NW, Airdrie, Alberta, T4B -- 100%

Why does the locator return a lower match score percentage for the "correct" street type and assign a higher match score percentage to an "incorrect" street type?

I can tell you from the outset, that the street data I'm using comes from a GIS vendor that has no intention of altering the way they designate their street types: Street types will always be 2 characters.


texts: How to get a native srs and declared srs of layer using gsconfig?

texts:
How we can get native SRS and declared SRS of existed layer in GeoServer?

I'm getting up to layer name available in a particular store and work space.

The code as follows

from geoserver.catalog import Catalog
cat=Catalog("http:// localhost:8080/geoserver/rest")
myResource = cat.get_resources(store='storeName', workspace='workspaceName')
layers = cat.get_layers(myResource)
for rss in myResorce:
print rss.name


texts: Enable feature access editing in Arcgis Server

texts: I had published a WFS/Feature Access service via ArcMap. But I can't edit layers (change geometry or data fields). Such settings in properties window are disabled (see an image). How to enable editing these layers?


texts: How to convert Latitude and Longitude to SPC

texts:
I am working with a set of data that is in Latitude and Longitude. I need to convert this data to State Plane Coordinates. Does anyone know a formula to be able to do this or a Java library that can do this?

I have been looking online have found this from NOAA but I am needing to implement this in Java for a piece of software and would need something with an API or just a formula.


texts: Can i use GPS Garmin GPSmap 62s in GIS projects with high accuracy?

texts: i have GPS Garmin GPSmap 62s, is this can be dependable in GIS projects with high accuracy, and what is the "RMS", the minimum error rate in meter, and can i linking it with arcpad to draw road or buildings with high accuracy.


texts: Converting CDF file to NetCDF file to use in ArcMap?

texts:
This is a GIS question as well as a Remote Sensing question.
I have downloaded CDF files from the ESA from MyEarth online and want to use them in ArcMap. In ArcMap you can Make NetCDF a Raster Layer. Is there some way I can convert CDF to NetCDF? Or do I need to? Should a CDF file work here?


texts: How determine optimal epsilon value in meters for DBSCAN by plotting KNN elbow

texts:
Before doing DBSCAN I need to find optimal epsilon value, all the points are geographical coordinates, I need the epsilon value in meters before convert it to radians to apply DBSCAN using haversine metrics
from sklearn.neighbors import NearestNeighbors
neigh = NearestNeighbors(algorithm='ball_tree', metric='haversine', n_neighbors=4)
nbrs = neigh.fit(np.deg2rad(firms[['y', 'x']]))
distances, indices = nbrs.kneighbors(np.deg2rad(firms[['y', 'x']]))

AND THEN
# Plotting K-distance Graph
distances = np.sort(distances, axis=0)
distances = distances[:,1]
plt.figure(figsize=(20,10))
plt.plot(distances)
plt.title('K-distance Graph',fontsize=20)
plt.xlabel('Data Points sorted by distance',fontsize=14)
plt.ylabel('Epsilon',fontsize=14)
plt.show()

and the graph output is this, but I need the epsilon value in meters.


texts:
I'm having an issue with ASTER data.

I have two mosaics; one for 2008 and one for 2012. Both are composed of 4 ASTER scenes and the scenes are all from the same two month period (i.e. no seasonal differences).

I have converted the values from DN values to Radiance (CORRECTION: reflectance value was calculated after radiance) and calculated for NDVI but the end result is absolutely meaningless (i.e. the NDVI range is between 1 and -1 but the comparison between two years can not detect change).

I think the reason for this is because the 2012 image looks very 'over-exposed' (i.e. there's little spectral difference between urban and pasture).

I think I might need to correct for atmospheric 'noise' but I had thought that AST14OTH images were already correct and that it wouldn't be necessary for ASTER14OTH images.

If anyone can shed light on why the images from the two years are so different and how to correct them - I would love to hear your advice.

Thanks!!!


texts:
I'm using ARC View v10 and have a requirement to buffer 1500+ building polygons by 2m. When buffered an overlap is created as shown highlighted in green in the image below.

Where as I would normally use a dissolve or merge to leave only the outline, I need to preserve the identity of each individual polygon after removal of the highlighted areas.

Can anyone out there demonstrate how this could be achieved?


texts: I basically use REST URL to display feature layer or ArcGISDynamic Service layer and it is commonly promoted by ESRI by their api example. I was wondering since the ArcGIS server also provide SOAP URL along with the REST URL so is it possible to use SOAP URL to configure feature layer or ArcGISDynamic Service  using javascript. So I can put SOAP URL to good use also SOAP protocols obviously has some good advantages.


texts: How do I convert tile pixels to longitude and latitude for web mercator?

texts:
I need a function to convert from tile pixels (x, y) from the Web Mercator projection (aka EPSG::3857 or SRID 900913), to longitude and latitude.

I'd like the function to be written in pseudo-code or a language like C.  I'm writing a C++ program.

I apologize if this is a duplicate question but I'm having trouble finding the algorithm.

I have found an algorithm that converts longitude and latitude to tile pixels.  But I'm not confident that I can rewrite it to do the opposite.


texts: Missed "processing' panel in QGIS 2.12.2 Lyon

texts: I have missing 'processing' panel in my QGIS Lyon version. It happened once some time ago, so I re-installed program then. And now I have the same situation -  I've cleaned my QGIS registry key and nothing changed. Anyone have any idea why? And how I can bring it back without re-installing?


texts: Identifying clusters in classified imagery data

texts:
I am relatively new to working with satellite imagery at 1 m resolution and am struggling with a methodology for "grouping" the imagery.

Using ArcGIS, I have performed an ISO classification on some aerial imagery and am trying to categorize different vegetation types. I am "see" the delineations between vegetation types, but they are based on multiple categories from the ISO classification. Combining the classes does not help. I have also worked with segmentation mean shift, but this yielded less than satisfactory results.

I am wondering if these is some methodology that allows for the reclassification of a classified raster that takes into account the class of neighboring cells. I know this is a relatively open-ended question, but not sure how to clear it up.

I think that the solution would be akin to classifying an impressionist painting.


texts:
I have a raster which is split into 30 tiles and I want to dump them all into PNGs.

I've read the documentation http://postgis.net/docs/RT_ST_AsPNG.html but still don't know how to dump it well.

If i use ST_Histogram, i get numbers, think the raster is well inserted,

SELECT band, (stats).*
FROM (SELECT rid, band, ST_Histogram(rast, band) As stats
      FROM "corine-g100_06-v15" CROSS JOIN generate_series(1,3) As band
      WHERE rid=2) As foo;


texts: Using arcpy.da.Walk to go through many file geodatabases

texts:
I have several file geodatabases I've downloaded from a content provider.  A subset of them, which I'm filtering by a couple of strings, have 1 feature dataset within (containing many feature classes) that I'd like to list and copy - at the feature dataset level - into one final geodatabase.

I am running into problems with getting the feature datasets isolated and listed.

Here is my code:

import os, zipfile, arcpy

pathPrefix = "C:\\Incoming"

newPath = '/FINAL'
try:
    os.makedirs(pathPrefix + newPath)
except:
    pass

for f in os.listdir(pathPrefix):
    <DO UNARCHIVING/EXTRACTION HERE>

Here's the part I'm having issues with:

fileFilter = ['_DA_10', '_DB_10']

for path, subdirs, files in os.walk(pathPrefix + newPath):
    for name in subdirs:
        if name.endswith(".gdb"):
            if any(x in name for x in fileFilter):
                print(name) ## This gives me the correct FGDBs I'm looking for
                workspace = os.path.join(pathPrefix, name)
                walk = arcpy.da.Walk(workspace, datatype = "FeatureDataset")
                for dirpath, dirnames, filenames, in walk:
                    for filename in filenames:
                        print filename

It runs, but gives me no feature datasets listed.


texts: How can i clip this TIN and form like the polygon.


texts: Conversion from N,E to Lat/Long And Lat/Long to N,E will always be same?

texts:
I have an excel tool for the conversion of following lat/long:

Lat: 23° 11’ 24.67503’’ S
Long: 118° 46’ 6.15742’’ E , which gives me following N,E:

Northing: 7,434,317.43

Easting: 680,991.69

But I am not clear that why conversion of the N,E back to Lat/Long is not giving me same results i.e. Lat: 23° 11’ 24.67503’’ S

Long: 118° 46’ 6.15742’’ E

Can anyone help me please?


texts: Trying To Identify The System That Coordinates Belong To

texts: I have architectural drawings of a bus terminal in Rochester, NY done in 1974 that lists the coordinates of the corners of each building. But I don't recognize the system of these coordinates. For example, instead of the usual "North" and "West," the first number is identified as "South," the second as "East." A full example is 49,876.06 South, 41,834.17 East. I need to convert these to lat/longs (degrees, mins, seconds) to relate to coordinates reported by the on-board GPS system of buses moving about the terminal. I thought these might be State Plane Coordinates, but I've found those have a different format. Can anyone identify the system these coordinates belong to?


texts: How to Predict demand on a bus route, considering variaton in kms and price?

texts:
If there is a route (line) modification, I mean change on distance; and also a change on the price for the user. How can I calculate the demand on the next year, considering seasonnality. What would be the best approach to solve this?

Multiple regression?
Time Series?
Spatio-Temporal predictions?

I want to obtain the result month by month


texts: How to fit census tracts within Voronoi Polygons?

texts:
I have a map with Voronoi polygons and I need to join census data  in it. The problem is that the geometry of voronois and census tracts (CT) do not match.

My goal is is to place census tracts (or only part of them, in most cases) within the corresponding voronois.

I'm having trouble thinking about how to calculate or weight the census data, since the census tracts might be contained in more than one voronoi polygon at the same time. Say, for example, how can I calculate the average population for any area of influence?


texts: Running PyQGIS script when opening any project in QGIS?

texts:
This question is a follow-up of this topic:
Setting Alpha Slider and Feature Count ON by default using QGIS?

I am trying to run this script for each new project I open.

But I am failing of making it active.

From my knowledge of QGIS, a script called startup.py should be created.

In QGIS 3, this script is there

[Linux] /$HOME/.local/share/QGIS/QGIS3/profiles/default/python

[Windows] C:\Users\You\AppData\Roaming\QGIS\QGIS3\profiles\default\python\startup.py

(You can check this directory under : Settings > User Profile > Open Active Profile Directory)

My startup.py script is :

from qgis.core import QgsMapLayer
from qgis.core import QgsProject
from qgis.utils import iface
import enableFeatCountAndAlphaSlider    #Activation sub

QgsProject.instance().legendLayersAdded.connect(enableFeatCountAndAlphaSlider)

The routine enableFeatCountAndAlphaSlider is stored in the same directory as startup.py and as the following content:

def enableFeatCountAndAlphaSlider(layers):
# By Joseph (https://gis.stackexchange.com/questions/329919/qgis-3-4-setting-alpha-slider-and-feature-count-on-by-default/329984#329984)
    root = QgsProject.instance().layerTreeRoot()
    layer = layers[0]
    # Enable feature count for vector-type layers
    if layer.type() == QgsMapLayer.VectorLayer:
        myLayerNode = root.findLayer(layer.id())
        myLayerNode.setCustomProperty("showFeatureCount", True)
    # Enable transparency slider
    if layer.customProperty("embeddedWidgets/count") != 1 or layer.customProperty("embeddedWidgets/0/id") != u'transparency':
        layer.setCustomProperty("embeddedWidgets/count", 1)
        layer.setCustomProperty("embeddedWidgets/0/id", "transparency")
    # Refresh legend symbology
    iface.layerTreeView().refreshLayerSymbology(layer.id())

# Connect "legendLayersAdded" event to "enableFeatCountAndAlphaSlider" function
QgsProject.instance().legendLayersAdded.connect(enableFeatCountAndAlphaSlider)

I have got no message from qgis windows.
When running manually startup.py from the console I have got the following error message:

  QgsProject.instance().legendLayersAdded.connect(enableFeatCountAndAlphaSlider)
  NameError: name 'QgsProject' is not defined

What am I missing and how do I correct it?


texts: Geoportal resource URL to Geoserver does not show Preview link in search result

texts:
I have a Geoportal catalog entry with a resource URL like:

  http://data.ilri.org:8080/geoserver/nurc/Arc_Sample/wms?request=getcapabilities

However the preview link is not shown in the search result.

I thought that maybe this was related to the resourceLinkBuilder.preview.filter variable in gpt.xml but even if I set the filter to empty (No filters at all)

<parameter key="resourceLinkBuilder.preview.filter" value=""/>

The preview link still does not show.

Any ideas how to make it work?

Please note that if I call the Geoportal Preview Window with my URL the Preview works fine.

  http://data.ilri.org/geoportal/catalog/livedata/preview.page?url=http://data.ilri.org:8080/geoserver/nurc/Arc_Sample/wms?request=getcapabilities

Many thanks,

Carlos


texts: Save the output in temporary memory under QGIS Plugin?

texts:
I am trying to make a simple plugin, from the input as a vector layer then processed with a Field Calculator function (But I want the ouput from this process to be saved in temporary memory or in other words the output to enter the next step function as an input) after that I want to do a Rasterized function.

I have tried to create the code :

def setVariables(self):

    self.rooftop = self.getRooftopLayer()
    self.rain = self.getRainfallLayer()
    self.potvol = self.dlg.le_potvol.text()

def Rasterize (self):

    alg_params = {
        'FIELD_LENGTH': 10,
        'FIELD_NAME': 'area',
        'FIELD_PRECISION': 3,
        'FIELD_TYPE': 1,
        'FORMULA': ' $area ',
        'INPUT': self.rooftop,
        'NEW_FIELD': True,
        'OUTPUT': QStringLiteral( "TEMPORARY_OUTPUT" )
    }

    outputs ['FieldCalculator'] = processing.run('qgis:fieldcalculator', alg_params, feedback=feedback)

    feedback.setCurrentStep(1)
    if feedback.isCanceled():
        return {}

    alg_params = {
        'BURN': 0,
        'DATA_TYPE': 5,
        'EXTENT': outputs['FieldCalculator']['OUTPUT'],
        'FIELD': 'area',
        'HEIGHT': 30,
        'INIT': None,
        'INPUT': outputs['FieldCalculator']['OUTPUT'],
        'INVERT': False,
        'NODATA': 0,
        'OPTIONS': '',
        'UNITS': 0,
        'WIDTH': 30,
        'OUTPUT': self.potvol
    }

    outputs['RasterizeVectorToRaster'] = processing.run('gdal:rasterize', alg_params, feedback=feedback)

def addLayers(self):

    self.iface.addRasterLayer(self.potvol, str.split(os.path.basename(self.potvol), ".")[0])

What statement I should add?


texts: Using ArcObjects to choose GeoTransformation?

texts:
I’m currently building a Desktop add-in tool with ArcObjects that:

Asks a user to select a feature class
Reprojects the feature class to Web Mercator
Executes some geoprocessing

The initial coordinate system of the feature class could be one of many different geographic or projected systems.  As a result I need to also have the user select a GeoTransformation if necessary.  Obviously, I could present the user with the huge list of transformation provided in the enumerations of esriSRGeoTransformationType, esriSRGeoTransformation2Type, esriSRGeoTransformation3Type.   But that would be a huge list.  What I’d like to do is narrow that list based on the input and output Coordinate Systems – but I have not been able to figure out how to do that narrowing.

Anyone have experience doing this? I know there must be some way to do, because the Project Tool UI does exactly this narrowing operation.  But I can’t find the method, despite an exhaustive internet search.


texts:
I need to remove Z dimension on spatial column (ST_GEOMETRY) with Esri on Oracle using the PL/SQL packages
I need to do the following, but using the ESRI PL/SQL packages on Oracle.
I have used a couple of the packages such as SDE.ST_GEOMETRY_OPERATORS, however I am unsure if doing the below (removing Z) is possible in PL/SQL.
There are directions of how to do this in ArcMap:
https://support.esri.com/en/technical-article/000010389

Browse to ArcToolbox > Conversion Tools > To Geodatabase.
Open the Feature Class to Geodatabase (multiple) tool.
Add all the feature classes into the Input Feature Class parameter.
Select an Output Geodatabase.
Click the Environments button at the bottom of the tool dialog box.
Expand the General Settings.
For the parameter, Output has Z Values, change the value to Disabled.
For the parameter, Output has M Values, change the value to Disabled.
Click OK in the Environments dialog box.
Click OK to execute the geoprocessing tool.

How can I do the equivalent in PL/SQL?


texts: Seeking border of Canada as GeoJSON from "official" source?

texts:
I'm new to GIS.  I need to get the border of Canada as a GeoJSON file.

I googled and landed at this repo https://github.com/johan/world.geo.json which as a representation that would suit my immediate needs, but it really feels like I should be able to go to "The Canadian Department of Interior" or something to get something a little "official".

Is there a canonical resource for this type of thing for Canada?


texts: Using directory names starting with number in ArcPy?

texts:
I have a script to export DDP extents to tif and when I change

arcpy.mapping.ExportToTIFF(mxd,"tests\BA\BA10k_" + str(ddp.pageRow.PageName) + ".tif",df,9306,9306*ar,1800,True)

to

arcpy.mapping.ExportToTIFF(mxd,"tests\BA\10k\BA10k_" + str(ddp.pageRow.PageName) + ".tif",df,9306,9306*ar,1800,True)

--added "10k" to the directory I get

  Exporting page 1 of 760 Traceback (most recent call last):   File
  "P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Creation\DDP_GeoTiff_production_BA.py",
  line 10, in
      arcpy.mapping.ExportToTIFF(mxd,"tests\BA\10k\BA10k_" + str(ddp.pageRow.PageName) + ".tif",df,9306,9306*ar,1800,True)   File
  "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\utils.py", line
  181, in fn_
      return fn(*args, **kw)   File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\mapping.py", line 1373, in
  ExportToTIFF
      layout.exportToTIFF(*args) AttributeError: DataFrameObject: Error in executing ExportToTIFF

Any ideas on why?


texts: Map material from China not allowed to leave that country

texts:
I just read an article about a company in Germany which tests autonomous driving cars in different countries.

They state that they analyze all of the generated data in Germany - with one exception: China.

  Because of legal regulations, no map material is allowed to leave the country.

Googling a few seconds for maps of china shows several results.

Can anyone explain this discrepancy?


texts: PYQGIS 3: Save custom layout in .qgs project

texts:
I iteratively create .qgs projects, load different layers in them, create a layout and export the maps. The .qgs projects are created but when I open them, my custom layout is not there! Therefore it was not saved within the project.
What did I miss?

Here is my piece of code:

processing_name = ["qgs_t1", "qgs_t2", "qgs_t3"]
anriss_list = ["C:\\MyData\\A300.shp", "C:\\MyData\\B300.shp", "C:\\MyData\\C300.shp"]
anriss_name_list = ["A300","B300","C300"]

for i in range(len(3)):
    # create project
    qgis_project_name = "C:\\ResultFolder\\" + processing_name[i] + ".qgs"

    project = QgsProject.instance()
    project.setCrs(QgsCoordinateReferenceSystem(21781))

    # add vector layer
    iface.mainWindow().blockSignals(True)
    shplayer = iface.addVectorLayer(anriss_list[i], anriss_name_list[i], "ogr")
    if not shplayer.isValid():
       print("Anriss .shp failed to load!")
    shplayer.setCrs(QgsCoordinateReferenceSystem(21781))

    project.addMapLayer(shplayer)

    #-- Create custom caneva
    manager = project.layoutManager()
    layout_name = 'my_custom_layout'
    layout_list = manager.printLayouts()

    layout = QgsPrintLayout(project)
    layout.setName(layout_name)
    layout.initializeDefaults()
    manager.addLayout(layout)

    #--- create map item

    map1 = QgsLayoutItemMap(layout) # create Layout item
    map1.setRect(QRectF(10, 10, 210, 90)) # set the map item rectangle
    map1.setPos(10,10) # map item upper left corner position
    map1.setFrameEnabled(True) # to have a map black border

    ms = QgsMapSettings()
    ms.setLayers([shplayer]) # set layer to be mapped
    map1_extent = QgsRectangle(shplayer.extent())
    map1_extent.scale(1.0)
    ms.setExtent(map1_extent)
    map1.setExtent(map1_extent)

    layout.addLayoutItem(map1)
    map1.attemptMove(QgsLayoutPoint(10,10))
    map1.attemptResize(QgsLayoutSize(210,90))

    # export map in .pdf
    exporter = QgsLayoutExporter(layout)
    exporter.exportToPdf("C:\\mapstest.pdf", QgsLayoutExporter.PdfExportSettings())

    project.write(qgis_project_name[i])
    project.clear()


texts: Creating scatterplot with icons (svg or png) instead of points?

texts:
I'm just starting with R and I wanted to know if it's possible to create a scatterplot with custom icons (svg or png) instead of points.

I started to import a shapefile, and then I have a set of points (latitude and longitude) that represents the location of the icons.

rm(list=ls())
library(ggplot2)
library(gstat)
library(sp)
library(maptools)
library(rgdal)

#plot
par(mar = c(2, 2, 2, 2))
 root <- "C:/"
Al= readShapeLines("/DZA_adm0.shp")
 require(rgeos)
simple.output <- gSimplify(Al, topologyPreserve = TRUE, tol = 3)
 plot(simple.output,col="black")
(tab <- read.delim(file = paste(root, "/ville.txt", sep = "")))
dat <- SpatialPoints(coords = tab[, 2:3])
(dat <- SpatialPointsDataFrame(coords = tab[, 2:3], data = tab[,
                                                               + c(1, 4)]))
# Rajout des capitales
plot(dat, pch =20, add = T,col="red")

# Extraction des coordonnees
lon <- dat@coords [, "lon"]
lat <- dat@coords [, "lat"]
# Extraction du nom des villes
labs <- as.character(dat@data[, "ville"])


texts: Removing commas from ArcGIS Online legend?

texts:
Right now I am working on a map using zip codes and have gotten the thousands comma out of the attribute table, but that is not changing the way it is showing up in the legend.

Does anyone out there know how to remove the comma from the legend?


texts: CartoCSS different background colors of Map element

texts:
I tried to change the background color of a map in zoomlevels greater than 11, but

Map{
  background-color:@water;
  [zoom>11]{background-color:@water*1.05;}
}

had no effect.
Is there an other way in CartoCSS to change the background-color for different zoom levels?


texts: Defining KML Name element when using QgsVectorFileWriter

texts:
I am trying to use QgsVectorFileWriter to write a kml file from a Shape file.
I would like to define the KML Name (called NameField in the QGIS save-as GUI) to correspond with a field in my shapefile, which is called "NAME".
The solution given in How to define KML Name element when using QgsVectorFileWriter? doesn't work for me.
Below is my code
from qgis.core import QgsVectorFileWriter, QgsVectorLayer

data_source = r"C:\Users\....\PT.shp"
layer = QgsVectorLayer(data_source, "layer_name", "ogr")
output_layer = r"C:\Users\....\PT_test"
QgsVectorFileWriter.writeAsVectorFormat(layer, output_layer, "utf-8",layer.crs(),"KML",datasourceOptions = "NAME")

print ("Layer to KML Conversion Complete")
print (output_layer + " is now a .kml")


texts: Postgis shapefile loader 2.1 error: dbf file can not be opened

texts:
The .shp, .dbf and the rest of files are fine, I open them in QGIS, AGOL with no problem. In fact, I've tried 3 different layers, and get always the same error.

Anyone can please help, as I really need to migrate my .shp's to postgis, and have no idea of why it's not able to open .dbf?

Postgres v9.3 x86
Postgis shapefile loader v2.1

Thanks a lot in advance!

PS I've tried changing the default character encoding from the default 'UTF-8' to 'ISO 8859-14' but same error.


texts: Why is Loop in Python Script Tool Validation only Looping Once?

texts:
My Python Script tool (ArcGIS Desktop 10.2) has two getParameterAsText inputs.

The first parameter is populated by the initializeParameters(self) module in the validation script which selects all CSVs in a particular folder.  Obviously set to 'Input' and 'Multivalue = yes'  It works fine.

The second parameter is populated by the updateParameters(self) module in the validation script by:

1) taking one CSV at a time (from parameter 1), and according to the information in its name...

2) pulling the proper path + file name from a "settings" (dbf) table with which to populate the second parameter, which is 'Output' and obviously set to "Multivalue = yes'.

The Problem:

If I select multiple CSVs from Parameter 1 the "Target_FCs" list only has the first corresponding path + file name appended to it (again, via the "settings" (dbf) table as described in step 2, above).

There are no problems with the values from my "settingsTable" (e.g. "mypath").

Quite sure the issue is that the for row in cursor line is only working the first time, just don't know why and how to make it work each time through the for csv in CSVs.split(';') loop.

BTW - I can confirm that the for csv in CSVs.split(';') loop is working correctly.  If I choose 3 CSVs from Parameter 1, this loop is traveled 3 times, but the for row in cursor is only traveled once.

Validation Code Snippet:

import arcpy
import os

  def updateParameters(self):
    settingsTable = os.path.join(os.path.split(os.path.realpath(__file__))[0],"Settings_CSV_to_ASMAD_FC.dbf")
    cursor = arcpy.SearchCursor(settingsTable)
    Target_FCs = []

    CSVs = str(self.params[0].value)
    x = 0
    if self.params[0].altered:
      for csv in CSVs.split(';'):
        if csv.upper().find("SMALLSITES") != -1:
          datatype = "SmallSites"
        elif csv.find("BIGSITES") != -1:
          datatype = "BigSites"
        if csv.upper().find("INACTIVE") != -1:
          activity = "Inactive"
        elif csv.upper().find("ACTIVE") != -1:
          activity = "Active"

        for row in cursor:
          myname = row.getValue("Name")
          mypath = row.getValue("Location")
          if myname == datatype + "_" + activity:
            Target_FCs.append(mypath)

        self.params[1].value = Target FCs

    return

CSV contents are not relevant to this question, only their names. Some names could be: "BIGSITES_ACTIVE.csv", "SMALLSITES_INACTIVE.csv", etc...

Settings Table looks like this:

Name                    Location
BIGSITES_ACTIVE         Q:\Data\Sites.gdb\BigSites_Active
BIGSITES_INACTIVE       Q:\Data\Sites.gdb\BigSites_Inactive
SMALLSITES_ACTIVE       Q:\Data\Sites.gdb\SmallSites_Active
SMALLSITES_INACTIVE     Q:\Data\Sites.gdb\SmallSites_Inactive

It's really just a lookup table to determine where the new FC (that is to be created from the CSV) is to be stored.


texts: Calculating urban areas that need pharmacies using QGIS?

texts:
I'm very new in QGIS.

The task is to calculate the urban areas that has lack of pharmacies.
I have buildings points with apartments count and existing pharmacies points. Then I made two raster heatmaps:

Buildings heatmap weighed by apartments count (blue on the image)
Pharmacies heatmap (red on the image)

I tried raster calculator with formula "apartments/pharmacies". The result looks like thruth (the more greener area, the more lack of pharmacies), but very ragged and exactly in bounds of pharmacies heatmap (not included other buildings areas)

And I cant go further.

What should be my direction?


texts: ArcObjects 10.1 Misplaced graticule labels after map export to jpeg

texts:
I´m currently working on an ArcObject 10.1 Extension which 'produces a Map' and exports it´s Page Layout to 'jpeg' after completion. The code follows the snippet for 'Create JPEG (high-resolution) from ActiveView':

     public static System.Boolean CreateJPEGHiResolutionFromActiveView(ESRI.ArcGIS.Carto.IActiveView activeView, System.String pathFileName, short shortQuality, int intResolution)
    {
        pathFileName = pathFileName.Replace(".jpeg", ".jpg");

        //parameter check
        if (activeView == null || !(pathFileName.EndsWith(".jpg")))
        {
            return false;
        }
        ESRI.ArcGIS.Output.IExport export = new ESRI.ArcGIS.Output.ExportJPEGClass();
        export.ExportFileName = pathFileName;

        IExportJPEG pExportJPEG = export as IExportJPEG;
        pExportJPEG.Quality = shortQuality;

        // Because we are exporting to a resolution that differs from screen
        // resolution, we should assign the two values to variables for use
        // in our sizing calculations
        System.Int32 screenResolution = 96;
        System.Int32 outputResolution = intResolution;

        export.Resolution = outputResolution;

        ESRI.ArcGIS.esriSystem.tagRECT exportRECT; // This is a structure
        exportRECT.left = 0;
        exportRECT.top = 0;
        exportRECT.right = activeView.ExportFrame.right * (outputResolution / screenResolution);
        exportRECT.bottom = activeView.ExportFrame.bottom * (outputResolution / screenResolution);

        // Set up the PixelBounds envelope to match the exportRECT
        ESRI.ArcGIS.Geometry.IEnvelope envelope = new ESRI.ArcGIS.Geometry.EnvelopeClass();
        envelope.PutCoords(exportRECT.left, exportRECT.top, exportRECT.right, exportRECT.bottom);
        export.PixelBounds = envelope;

        System.Int32 hDC = export.StartExporting();

        activeView.Output(hDC, (System.Int16)export.Resolution, ref exportRECT, null, null); // Explicit Cast and 'ref' keyword needed
        export.FinishExporting();
        export.Cleanup();

        return true;
    } // end of CreateJPEGHiResolutionFromActiveView

After export it occurs from time to time that the graticule labels are displaced, whereas they are at the correct position in the Layout of ArcMap.

Does anybody has a hint for me what could be the cradle of this mess and how to work around this Problem?