texts: Handling MODIS NoData cells with gdalwarp

texts:
I'm using the gdalwarp tool to convert MODIS .hdf into GeoTiff files. I want to set the value of the nodata cells so that when I import the .tif image into Grass there is no black edge. I'm using the following command:

gdalwarp -of GTiff -t_srs EPSG:4326 -dstnodata 0 <hdf> out.tif

This is working fine but I'm concerned that by setting -dstnodata to 0 I may be losing the valid cells that happen to have a value of 0. Is this actually what is happening, or does the tool work in a different way?

Thanks!


texts: getGetFeatureInfo with time parameter returns value of old layer

texts:
I have a problem when I try to get a value from my group of layers.

I have some NetCDFs in geoserver that update every 6 hours.

In my code, I get the data from them with getGetFeatureInfoUrl. I use the parameter time to specify which value I want.

http://server/geoserver/trafikoa/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=trafikoa%3Aval_met&tiled=false&STYLES=&LAYERS=trafikoa%3Aval_met&tilesOrigin=-11.997193701565266%2C35.905838049948215&INFO_FORMAT=text%2Fjavascript&FEATURE_COUNT=50&X=117&Y=65&WIDTH=256&HEIGHT=256&SRS=EPSG%3A3857&BBOX=-303302.12823557854%2C5283327.395071384%2C-298410.1584253273%2C5288219.364881635&time=2017-02-17T00:00:00.000Z

The first time I get the value, it returns OK, but when I update the NetCDF, it returns the value from 06:00:00 and not 00:00:00.

I think, somehow it's cached or something because if I remove the layer from geoserver and create again, it works fine.

How do I have to set up Geoserver?


texts: Projecting geometry features upside down using MapServer

texts:
I have shapefiles containing polygons and lines that are used in a mapfile for projecting them through Mapserver.

Is there a way in changing the projection of the layers by 180 degrees (i.e. upside down)?


texts:
I have an openlayer map on the website.
I need to draw lines, points and polygons, but after it is draw, the user should be able to change the position of the geometry.
In my map it is already done the part of the drawing geometry type, the part with modifying the geometry type is missing.
I've checked this example: https://openlayers.org/en/latest/examples/draw-and-modify-features.html and this one: http://tsauerwein.github.io/ol3/animation-flights/examples/draw-and-modify-features.html.
It is explained very well, but since I have no experience with open layer, it is very hard to use this example in my case. I always get errors, or the map is not loaded at all.

I'm using this code to draw a point, line or polygon:

function createDrawForPoint(type) {
draw = new ol.interaction.Draw({
    source: source,
    type: type
});

draw.on(
    'drawstart',
    function (evt) {
        sketch = evt.feature;
    },
    this
);
draw.on(
    'drawend',
    function (evt) {
        var message = '';
        var title = '';

        geom = sketch.getGeometry();
        switch (type) {
            case 'Point':
                message = ol.coordinate.format(ol.proj.transform(geom.getCoordinates(), 'EPSG:900913', 'EPSG:31700'), '[{x}, {y}]', 0) + ol.coordinate.format(ol.proj.transform(geom.getCoordinates(), 'EPSG:900913', 'EPSG:4326'), '[{x}, {y}]', 4);
                title = "POINT [Stereo 70 / WGS 84]";
                break;
            case 'LineString':
                message = geom.getCoordinates();
                title = "LINESTRING";
                break;
            case 'Polygon':
                message = geom.getCoordinates();
                title = "POLYGON";
                break;
        }
        //alert(message);
        $('#geometry-dialog').find('.modal-title').html(title);
        $('#geometry-dialog').find('.modal-body > p').html(message);
        $('#geometry-dialog').modal('show');
    },
    this
);

return draw;

}

I also have these lines of code where the function createDrawForPoint is invoked:

    'addpoint': {
    'obj': createDrawForPoint('Point'),
    'on': function () {
        map.addInteraction(this.obj);
    },
    'off': function (event) {
        var i,
            ov = map.getOverlays();

        sketch = null;
        for (i = 0; i < ov.length; i++) {
            map.removeOverlay(ov[i]);
        }
        for (i = 0; i < mapOverlays.length; i++) {
            map.removeOverlay(mapOverlays[i]);
        }
        mapOverlays = [];
        source.clear();
        map.removeInteraction(this.obj);
    }
},
'addline': {
    'obj': createDrawForPoint('LineString'),
    'on': function () {
        map.addInteraction(this.obj);
    },
    'off': function (event) {
        var i,
            ov = map.getOverlays();

        sketch = null;
        for (i = 0; i <= ov.length; i++) {
            map.removeOverlay(ov[i]);
        }
        source.clear();
        map.removeInteraction(this.obj);
    }
},
'addpoly': {
    'obj': createDrawForPoint('Polygon'),
    'on': function () {
        map.addInteraction(this.obj);
    },
    'off': function (event) {
        var i,
            ov = map.getOverlays();

        sketch = null;
        for (i = 0; i < ov.length; i++) {
            map.removeOverlay(ov[i]);
        }
        source.clear();
        map.removeInteraction(this.obj);
    }
},

Can you give me some hints how can I solve this problem?


texts: Using area() or $area function in QGIS, when the data (polygon) is stored in UTM?

texts:
I have a question regarding the area() and $area function in QGIS and St_Area in PostGIS.

I am trying to calculate the exact area of Polygons stored in a Shapefile in UTM 25832 and do not know which functions are the correct ones to use. I know that the difference between area() and $area in QGIS 3.10 is that the first is calculated planar and the second with ellipsoid.

I have also tried to calculate the area with ST_Area in PostGIS and get the same result as in QGIS with the area function or with ArcMap.

When my data is projected as UTM 32N (EPSG:25832) isn’t it already planar and can be calculated with the simple area() function instead of $area? Or do I have to use the $area function even if it is in UTM32?


texts: Downloading buildings data from OSM with a polygon (shapefile) as the bounding box

texts:
I am struggling to complete this task. I am trying OSMnx which can be used to download data from OSM above however I am getting an error when using its from_polygon feature when trying to download data. Also, I am not sure if this data will include Buildings data.
I load my shapefile into geopandas and can then view and interact with it
This is the code
Building_data = ox.graph_from_polygon(my_shapefile, network_type='all')
ox.plot_graph(Building_data)

However, I am getting this error
ValueError: The truth value of a Series is ambiguous. Use a.empty, a.bool(), a.item(), a.any() or a.all().

EDIT: I managed to convert the shapefile into a shapely polygon and run the OSMnx script again,
from shapely.geometry import shape
import fiona
c = fiona.open('E:/shapefile.shp')
pol = c.next()
geom = shape(pol['geometry'])

import osmnx as ox
data = ox.graph_from_polygon(geom, network_type='all')
ox.plot_graph(data)

however, now I get this CRS error:
Invalid projection: +proj=utm +zone=80957 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +type=crs: (Internal Proj Error: proj_create: Error -35 (invalid UTM zone number)).

It seems a shapely polygon has no CRS or the wrong CRS.


texts: Clipping a raster by a polygon layer: Error "Cutline is not valid after transformation"

texts:
I have a raster layer (3 cm x 3 cm) and a vector layer with polygons and when I try to clip the raster by the polygon layer to have just the bare earth/ground regions left over, so I can interpolate them afterwards, I always get the error below:

Warning 1: Ring Self-intersection at or near point 571.37754947319627 322.188317745924
ERROR 1: Cutline is not valid after transformation

The polygon layer was derived from the raster layer.

I am working with QGIS.


texts: Incorrect Routing Times in ArcGIS Desktop 10.7

texts:
I'm trying to find an OD Cost Matrix for a sample of zip codes in LA, San Bernardino, and Clark (NV) Counties using TIGER Data found here: https://www2.census.gov/geo/tiger/TIGER2012/ROADS/ (the files I used are 06037, 06071, and 32003)-- but the times are 50%+ off from Google maps.
Specifically, I used the Integrate tool on all three, with a clustering tolerance of .001 decimal degrees, since the original shapefiles have gaps and so can't perform routing analysis. (Integrating by 100m made little difference).
I then used Project to change the shapefiles' coordinate systems to "USA_Contiguous_Equidistant_Conic" and outputted to feature classes in a feature dataset.
Now that the integrated and projected shapefiles are in a feature dataset, and because they're in projected coordinates, all the lengths of the roads are listed in meters. I created a new field that divides that length by 1609.34 to get the length in miles. I also assigned each road a speed in miles/second based on its road type in a new field. (e.g. a highway, as defined by the Census, has a speed of 70/3600 miles/sec). The VB code I used was:
speed = 25/3600
Select Case [RTTYP]
Case "C"
speed = 40/3600
Case "I"
speed = 70/3600
Case "M"
speed = 30/3600
Case "O"
speed = 30/3600
Case "S"
speed = 50/3600
Case "U"
speed = 55/3600
End Select

All this created the following attribute table for each feature class:

I then created a network dataset, with global turns and an attribute for "time" in seconds, which is just the length of each road in miles divided by its speed in miles/second:

I then used a few zipcodes from my own data to create a 15 x 15 ODCM using the Network Analyst tool box.
As I've mentioned, this is way off from Google Maps, and upon closer inspection it seems like everything is connected, but the program simply ignores the most optimal routes. For instance, ArcMap says to drive from Las Vegas to Los Angeles requires driving on unnamed dirt roads in the desert. The difference in the routes are apparent especially in the beginning of the route near Vegas:

I also tried creating a hierarchy where roads the Census calls primary roads are primary, secondary are secondary, and all the rest are tertiary. However, when I did this, ArcMap couldn't find any routes! I am also sure that the primary roads are connected; I created a topology using the three shapefiles with the rule "Must Not Have Dangles". It seems like most of those errors occur when roads actually end, and so there seem to be no discontinuities:

Earlier, when I projected and then integrated (still by .001 decimal degrees), the errors were about the same. Interestingly, when I try to integrate before projecting, ArcGIS often throws a topological error, noting that there's an "undetected intersection." Integrating by .01 decimal degrees only made the errors worse.
I also tried this method on the shapefile for Cook County, IL, and its errors compared to Google Maps were only about 3%. Finally, I did the same using the "primary and secondary" roads in Indiana and Illinois, and the average error was only 5%.
I have no idea what happened. Are there other data sources which don't require messing around with the Integrate tool (I tried using the OSRM API, but it seemed inaccurate compared to Google Maps)?


texts: Bug: QGIS freezes editing large polygon, does not save window or toolbar positions after crashes

texts:
Version: 3.14.15
Activity causing freezing: attempting to delete more than 20 vertices of a 61,000+, complexly-shaped polygon - hangs when deleting as few as 4 vertices
Annoyance: layer and browser windows and default taskbar positions always reset to defaults; in addition, most recent project disappears from main starting window


texts: Possible range for MSAVI and MTVI indices

texts:
What are possible value ranges for MSAVI2 and MTVI indices?
The article #1 says it should be in range of -1 to 1, but the answer to the question #2 says '[MSAVI2] is not bound to a theoretical -1 to 1 range'. Ranges for a bunch of vegetation indices are given on page #3, but not for MSAVI2/MTVI.
For example, if I get stats like these for MTVI:
MAX=18
MIN=-50
MEDIAN=0.06
MEAN=0.15

and
MAX=0.99
MIN=-20
MEDIAN=0.22
MEAN=0.3

for MSAVI2 - does that indicate I'm doing something wrong and/or that the results are messed up?
If I try to visualise those results with min-max set to -0.3 to +1, the image looks very much as expected for a vegetation index. NDVI gives result with all the data ranging approximately from -0.99 to +0.99 on the same scene.
Data is Sentinel-2 image, formulas are exactly as in the answer form the link #2:
msavi2 = (2 * nir + 1 - sqrt( (2 * nir + 1)^2 - 8 * (nir - red) )) / 2
mtvi = 1.5 * (1.2 * (nir - green) - 2.5 * (red - green)) / sqrt( (2 * nir + 1)^2 - (6 * nir - 5 * sqrt(red) - 0.5) )
Links:

https://wiki.landscapetoolbox.org/doku.php/remote_sensing_methods:modified_soil-adjusted_vegetation_index#discussion_comments

Soil Adjusted Vegetation Index Values (MSAVI2) are Not Within Expected Range

https://www.researchgate.net/post/What_is_value_range_for_different_Vegetation_Indices


texts: Python - read shapefile .dbf for encoding

texts:
I've written a python script that formats a .shp such that its fields structure matches my office's database view in order to ingest the data of said shapefile.  This original shapefile was encoded with the regional Hangul (Wansung) encoding, which required the codecs.decode('949') method to read in the Korean characters.  This is consistent with the ESRI table found here.

What I'd like to do is enhance the script so that it can read the .dbf header, find the encoding, which according to this and other sites resides at the byte offset of 29, so that our staff does not have to know the ESRI language code for every regionally-encoded shapefile we get.  It's a one-byte flag, and I figured '949' would be an integer, so python's struct.unpack method should do the trick, specifically this:

myBinaryFile = angusMcFile.read(30)
dbfHead = struct.unpack('<29xB', myBinaryfile)

The returned values are 'W' when I unpack guessing that this 'flag' is a 1-byte string, and '87' when I unpack guessing it's a 1-byte integer.  This is not the '949' or '0x79' that I'm looking for.  What am I doing wrong? I'll be happy to provide more code if necessary.


texts: Failure to Revoke ArcGIS Pro Advanced license using Python API

texts:
I am scripting the provisioning and de-provisioning of users for my place of employment. I have the script in place to license users but I am having trouble getting the revoke script to work. I can pull the necessary parts out and run them manually in a python terminal and it still fails to revoke the license. When I run it in the terminal the only message I see is "False". When I run it in the script it still fails but apparently not enough to throw an exception. The account that is running this script can provision licenses just fine I have even given the account full administrator access and it still fails. I can log in with the account and de-provision licenses manually but it fails when scripted.

The other issue that I run into is that it won't let me de-provision all of the entitlements at once using this line if pro_license.revoke(username=user.username, entitlements= '*'):. I receive this error message (assign() got an unexpected keyword argument 'supress_email')

This is the snippet from my script:

for user in new_users:
try:
    pro_license = gis.admin.license.get('ArcGIS Pro')
except:
    logger.warning("Cannot query ArcGIS Pro licenses")
    exit(1)

try:
    if pro_license.revoke(username=user.username, entitlements= 'desktopAdvN'):
        logging.info("Removed license from : " + user.username)
        print("Removed license from : " + user.username)
    else:
        logging.warning("Failed to remove license from: " + user.username)
        print("Failed")
except Exception as e:
    logging.warning("Failed to remove license from: " + user.username + "(" + str(e) + ")")`

Am I missing something obvious?


texts: Remove unwanted regions from map data QGIS

texts:
I downloaded Census Tract shape file data from the US Census for the state of Maryland. There are ~1400 census tracts. This is a vector layer in QGIS.

I have a separate CSV that has 84 census tracts. I only want these 84 tracts to show and remove the rest of the census tracts. I've tried filtering in the attributes table but  haven't found a way yet. Is there a way to match the list of 84 census tracts with the ~1400 census tracts shape file and remove the rest?


texts:
We are currently getting some new survey marks established in the main town of our shire, we have been told that there are different orders based on accuracy: 1st, 2nd and 3rd both for horizontal and vertical.   In the main town we already have a few good 1st order horizontal but very few 1st order vertical marks, and we where told that they would just shoot from these marks to make the new ones.

My question is:  How do they establish a 1st, 2nd or 3rd order mark in a town that doesn't already have a good mark to use; and how do you know it's right?

I know this might not be directly GIS related but to me surveying is important to GIS.


texts: Are there are shortcuts for exporting polygon vertex location as a table

texts:
Someone came to ask me how he can export the X and Y locations of vertices from his polygon shapefile to excel.

So I told him to convert it to points, add columns, calculate geometry, and then export the attribute table. This isn't too difficult but I found it was quite a long explanation to give for what seems a simple task.

Are there any shortcuts like a tool that will do this in one step? I had a tool through xtools pro and ET Geo-wizards but didn't find anything obvious.


texts: How to find the watershed for each outlet point and sum values within it?

texts:
I have the following files to work with:

Stream network of a basin;
Outlet points for each segment* of the network and
Points which represent water withdrawals along of the stream, which are not located in the stream, but near of segments.

*segment = a line with a start point and an outlet point.

My goal is to calculate the sum of withdrawals located upstream of each outlet points.

My first attempt was using model builder, iterating over each outlet point. For each of them, I delineated the watershed and calculated the sum using Zonal Statistics (summing withdrawals contained by the watershed). But it didn't work well because ArcGIS aren't calculating watershed properly. I do not know why...

I'm wondering if I can figure out how to select directly upstream watersheds each iteration, assuming I've already obtained a feature containing the watersheds for each segment.

NOTICE: I don't want to find the watershed for each segment. I want to find the watershed for each outlet point. That's more comprehensive.


texts: What is the annual data volume produced by the individual Sentinel satellite missions?

texts:
I am wondering how big the volume of the data produced by Sentinel-1, Sentinel-2 etc. in a given year. For Sentinel-1, it'd be most interesting to know the data volume of all SLC images, for Sentinel-2 the data volume of all granules.

Does anybody know a source where one can find these numbers?


texts: Why is loading specific WMS giving me CORS error, while others work

texts:
A client wants to dynamically add and remove WMS Layers to a web application from a governmental source.

The code to do this is a simple as it can get:

//layerObj has a property id, that contains the url as a string
var layer = new WMSLayer(layerObj.id);
//wait until the infos have been retrieved
layer.on("load", function(i) {
    var visibleLayers = layer.layerInfos.map(function(x) {
        return x.name;
    });
    layer.visibleLayers = visibleLayers;
    //update the store and add the layer to the coresponding item, so we can later on remove it form the map, without
    //actually querying the map
    var item = _store.query({ id: layerObj.id })[0];
    item["layer"] = layer;
    item["title"] = layer.title;
    _store.put(item, { overwrite: true });
    //add to map
    map.addLayer(layer);
    //add to list
    addToLayerList(layer);
});

I am using http://www.geodatenzentrum.de/geodaten/gdz_rahmen.gdz_div?gdz_spr=deu&gdz_akt_zeile=2&gdz_anz_zeile=5&gdz_user_id=0 as a source to test different WMS. The first table under "Beispiele" gives you an actual url to implement.

http://sg.geodatenzentrum.de/wms_dtk200 is working
http://sg.geodatenzentrum.de/wms_dtk1000 is working
http://sg.geodatenzentrum.de/wms_dgm1000 is working

but once I get to http://sg.geodatenzentrum.de/wms_webatlasde.light

I am receiving the following error:

XMLHttpRequest cannot load http://sg.geodatenzentrum.de/wms_webatlasde.light?SERVICE=WMS&REQUEST=GetCapabilities.
The 'Access-Control-Allow-Origin'
header contains multiple values 'http://localhost:8080, *', but only one is
allowed. Origin 'http://localhost:8080' is therefore not allowed access.

Why would this work for most of the WMS available, but fail for one or 2?


texts: how to Find n nearest lines for the given points from a collection of lines?

texts:
I have two tables, points and lines. For each point in point table, I want to find out nearest 3 points in the nearest 3 cables.

I almost developed a query, but it is best suitable for single point, but i need a query for all the points in the point table.

My query is

select * from (
with index_query as (
select l.data,p.data, ST_Distance(c.geom,p.geom) as distance,
ST_AsText(st_line_interpolate_point(l.geom,
st_line_locate_point(l.geom, p.geom))) as nearest_point_in_line,
from lines l,points p order by l.geom <#> p.geom limit 100)
select *,RANK() OVER (PARTITION BY pid order by distance ) AS r
from index_query ) as s where r < 4

LIMIT 100 is the limiting factor, which is limiting total output rows to 100 irrespective of number of points in Point table.

I need nearest 3 point_in_line for all points in my table. Total 3*no of rows in point table. How to achieve it?


texts: Putting random points only on selected parts of a line in QGIS?

texts: I am trying to create random points along a line, however I only want to do it in the selected area. Is there a way to do this?


texts:
I am following this tutorial https://docs.djangoproject.com/en/1.8/ref/contrib/gis/tutorial/ and did layer mapping. And going along with documentation, I tried the below code

from world import load
load.run()

And i got an error in shell which is attached below. Why is it so?
I could guess it is due to the NULL value in the attribute. I wanted to know anyone who followed the same documentation got the same?


texts: Adding nodes to pipe ends using ArcGIS Desktop?

texts:
I have two features, a pipe distribution network, and nodes. I need to have nodes at the end of each pipe. I have already found which pipes have no nodes at their ends.

How can I add nodes to the end of the pipes that have no nodes?

I am using ArcGIS Desktop 10.4.


texts: Is there a way to QGis read an SDE.ST_GEOMETRY field?

texts: this is the type and name that Qgis displays for me when I connect with a oracle table,it is not displayed as a geo column, there is a setting I'm missing? Can Qgis Read this type of data?


texts: geoExt2 application in a CMS (joomla, Wordpress)

texts:
I am working on geoExt2 application but I want to integrate it as "article" in a Joomla (or similar CMS). I don't want special interactions with the CMS just a link from the menu which will serve the geoExt2 gis application.

Have anyone tried this ?
Are there ways of doing this or CMS & geoEXT2 are mutually exclusive ?


texts: Layer error: ImageCollection.mosaic: Invalid type. (Earth Engine)

texts:
An error occurred in raster to vector conversion:

  Layer 1: Layer error: ImageCollection.mosaic: Invalid type.
  Expected type: Image.
  Actual type: FeatureCollection.

Why is the Map function trying to mosaic an image collection when I have a feature collection that is also a valid input to Map.

My code link is: https://code.earthengine.google.com/f14646de49f536d7ed4b4e0d5a9ccb21

  var elevation=SRTM;
  var sorted = L8.filter(ee.Filter.eq('WRS_PATH', 147))
  .filter(ee.Filter.eq('WRS_ROW', 38))
  .filter(ee.Filter.calendarRange(2013,2013,'year'))
 .filter(ee.Filter.calendarRange(7,10,'month'))
 .filterMetadata('CLOUD_COVER', 'less_than', 10)

var ndwi1 = sorted.map(
 function(img) {

   var ndwi = img.normalizedDifference(['B3','B5']);
   var slop = ee.Terrain.slope(elevation);
   var evi1 = ndwi.expression('(thr>0.25) && (slp<20)? 1 : 0' ,

      {  'thr' : ndwi,  'slp' : slop   });

   var evi = evi1.clip(table)

   var vectors = evi.addBands(evi)
     .reduceToVectors({
      geometry: table,
      crs: evi.projection(),
      scale: 30,
      geometryType: 'polygon',
      reducer: ee.Reducer.mean(),
      maxPixels: 1e8
                     });

 return vectors;
 })

 print(ndwi1)
Map.addLayer(ndwi1)


texts: Avoid creating white margins in tile cache in ArcGIS 10.3

texts: I'm relatively new in Geographic Information Systems and in ArcGIS software in general so I have no idea how to handle the following situation. I'm using a ArcMap software to create and publish a map service, which should be used as a tile cache for our software. I've created a tile scheme that suits our organization the best, set the correct projection, and level of details. The problem I'm having is that once the service is published and cache is created, there are white margins on tiles on the edges. This is a problem since we use this tiles in our software and the image is shown shifted to the side by those margins. The map is pretty simple, it's just one georeferenced image, nothing special for now. Is there any way to remove these margins from tile images?


texts: How to find the direction of an “L” shape in a image

texts:
I  need to find the moving direction of a vehicle by its extracted point cloud, and I have converted the point cloud to the following image.

As the target vehicle could be moving straight or turning and the image is sometimes clear and sometimes fuzzy, I find it's difficult to match the "L" shape using template matching.

I also try to use RANSAC to fit the linear, but it has two sides and RANSAC does not work well. What I need to do is using an oriented bounding box to represent the vehicle.

If I could have the yaw angle of the "L" shape, it's very easy to recover it to an oriented bounding box. So could anyone give me some suggestions?

PS: I have asked this question at stackoverflow, but nothing helps. So forgive me again post this question here.
PS: The function cv::minAreaRect could offer a basic result, but it sometimes fit the "L" shape in a wrong direction.


texts: Adding Edit with Pyscripter back to Windows context menu?

texts:
Using Windows 10 64 Bit, ArcGIS for Desktop Advanced, ArcGIS Pro 1.2

I use PyScripter pretty extensively, and after I installed ArcGIS Pro (which installs Python 3), I noticed that I can no longer right-click on a python file and choose "Edit with PyScripter" anymore, even though the setting in PyScripter is set properly (Tools ==> Options ==> IDE Options ==> Shell Integration ==> File Explorer Context Menu)

How do I get this function to work again?


texts: Defining NoData values in ArcPy Clip raster tool?

texts:
I'm clipping some rasters with the ArcPy Clip tool (from the raster toolset) and a polygon feature and I'm having some problems with nodata values.

The script I'm using is based on the code from Arcgis help:

  Usage:
  Clip_management (in_raster, rectangle, out_raster, {in_template_dataset}, {nodata_value} {NONE | ClippingGeometry})

##Clip Raster Dataset with feature geometry
arcpy.Clip_management("image.tif", "#", "clip.tif","feature.shp", "0", "ClippingGeometry")

If I use '0' as NoData value, the clipping is OK. But there are some pixels where a band has 0 values and I get this in the clipped pixels:

So, despite the clipping area it's fine, the raster doesn't look like the original because there's a band without its value.

I have been trying to use the ExtractbyMask tool but I've seen that some raster values are modified. It doesn't works for me because I need exactly the same values from the source.

Does somebody know if it's possible to define in ArcPy a Nodata value that considers the 3 bands at the same time?

Is it possible to define NoData value as '0,0,0' to solve this question?

I've been trying this but I'm not used to Python coding and maybe I'm not doing it properly.


texts:
I am trying to use the command las2dem in terminal (Ubuntu) to create a DEM with various las files but I only want to conserve the points classified as Ground. I'm using this command:
 las2dem -i file.las -o dem.png -keep_class 2

It returns the error: las2dem:command not found

I installed LAStools with this command: conda install -c conda-forge lastools

I have used other tools from LAStools like txt2las and las2las and I haven't had any problem before, I think there is not a problem with the installation but I'm not sure.

I can use the tool in qgis but I only can add one las file I need to incorporate to the DEM at least 20 las files.

Is there any way to use las2dem from the command line?


texts: Faster way to calculate slope between two points on DEM using ArcGIS for Desktop?

texts:
I'm a forester, and I do lots of work with LiDAR derived DEM's. I'm looking for a quick-and-easy tool to calculate the slope between two points. The application for this would be drawing a proposed road, where the slope of each road segment must be below a certain value. After quite a bit of searching, I'm unable to find a tool to accomplish this in ArcMap. Can anyone suggest an existing tool to accomplish this?

Edit: I don't need the slope along the surface of a DEM; I simply need a FAST way to calculate the straight-line slope between two points. Currently, the only way I'm able to do this is as follows:

1) Load DEM into ArcMap.
2) Create a new shapefile.
3) Create a line feature with a single segment.
4) Get the elevation values for the lines two vertices from the underlying DEM.
5) Use math to calculate the slope of the straight line between the two vertices.

I'm looking for a workflow similar to the following:

1) Click on "slope" tool.
2) Click on point A, then click on point B.
3) Dialogue box displays the straight-line slope between point A and point B.

I've been looking for an existing solution to accomplish this, and have been really surprised that I've come up empty handed. Can anyone recommend a good way to do this?


texts: Why are the tile coordinates different in Google maps and Google earth?

texts:
I have a .NET application that typically uses the standard Google Maps JavaScript API V2. Within this application I am using a GTileLayerOverlay to add my own custom dynamic tile layer.

  function addJourneyLineOverlay(assetId, startDate, finishDate) {

        var path = "http://<%= Request.Url.Host %>";
        if ("<%= Request.Url.Port %>" != 80) {
            path += ":<%= Request.Url.Port %>";
        }
        path += "/Private/polyLineTileLayer/Tile.aspx";

        var myCopyright = new GCopyrightCollection("(c) ");
        myCopyright.addCopyright(new GCopyright('Demo', new GLatLngBounds(new GLatLng(-20, 130), new GLatLng(-30, 155)), 0, '@2012 IMT Pty Ltd'));

        var journeyPathLayer = new GTileLayer(myCopyright);
        journeyPathLayer.getTileUrl = function (tile, zoom) {
            return path + '?zoom=' + zoom +
              '&x=' + tile.x +
                  '&y=' + tile.y +
                      '&assetId=' + assetId +
                          '&startDate=' + startDate +
                              '&finishDate=' + finishDate;
        };

        journeyPathLayer.isPng = function () {
            return true;
        };
        journeyPathLayer.getOpacity = function () {
            return 1.0;
        };

        journeyPathOverlay = new GTileLayerOverlay(journeyPathLayer);
        map.addOverlay(journeyPathOverlay);
    }

This overlay requests it tiles from an aspx page that gets some points from a database, determines if the points are in the tile and then returns an image with the points on it. I developed it based on the same pattern as Gheat.net. And I am also using the Gmap.Net library for its Mercator projection, for when I need to work out which points should be in which tiles.

This works really well when I am using the standard Google maps API. However, recently I have to integrate my application with our clients Google Earth Enterprise server. This means that I am now getting my map imagery and API from our clients server and not from Google proper. This has caused a problem with my overlay.

My overlay is not longer position where it should be relative to the map location. They are appearing a bit lower that where they should be.

Using firebug I examined the request for my tile layer and found out that the tile X and Y where different than what I was expecting.

At zoom level 12, using Google maps a particular LatLng location is in tile X-3392, Y-2312. But in Google Earth the same location at the same zoom level is in tile X-3392, Y-2305. that's a 7 tile difference in the Y axis.

What can account for this difference? Does Google Earth use a difference projection, if so is there a convenient .net library for it?


texts: How to fill "shapefile name" and "name of attribute"?

texts:
When I try to create a shapefile two of the boxes to fill in are "name of the attribute" and " shapefile name".

This may be very basic but I can't find and answer.  I want to fill them in correctly.  What are their purposes and what should one put in them?  Why are their two names for the file?


texts: Cannot execute gdalwarp in a Python notebook

texts:
I try to reproject a raster using gdalwarp.
To do so I use the subprocess package of Python3 :
command = [
    'gdalwarp',
    map_raster,
    map_raster_proj,
    '-t_srs', 'ESRI:54009'
]
print(' '.join(command))

kwargs = {
   'args' : command,
   'cwd' : os.path.expanduser('~'),
   'stdout' : subprocess.PIPE,
   'stderr' : subprocess.PIPE,
   'universal_newlines' : True
}

with subprocess.Popen(**kwargs) as p:
    for line in p.stdout:
        print(line)

unfortunately I get the following error :

ERROR 1: PROJ: proj_get_authorities_from_database: Cannot find proj.db
ERROR 1: Translating source or target SRS failed:
ESRI:54009

When I run the command directly from my terminal It works flowlessly.
I imagine that I have some environment parameter missing or pointing to the wrong direction but I don't know where to start.
Any idea where I should check differences between the two execution env?
EDIT
I run Jupyter notebook from a company cloud interface, I have 0 idea of the installation, I have 0 right to access it. It's running on linux AWS instances.
several test that I performed :
!gdalwarp --version

give me the same result in the terminal and the jupyter notebook : GDAL 3.0.4, released 2020/01/28
!printenv | grep PROJ_LIB

I have nothing in the terminal env and PROJ_LIB=/usr/local/lib/python3.6/dist-packages/rasterio/proj_data in the notebook.
If I run this command in another notebook I get 0 errors and the PROJ_LIB variable is not set. I think my install of rasterio and gdal are coliding is it possible ?


texts:
Are you aware of any software/tools which permit the automatized or semi-automatized creation of old-fashioned maps based on data one owns?
I love the design of old-fashioned maps, in particular the symbols used for buildings, mountains and forests. They transfer a lot of visual information, with often limited use of colour.
A related possibility is the Stamen Watercolor style for basemaps, which looks very nice, but does not visualize a lot of information.
Here are some examples:
http://cdn.supadupa.me/shop/14281/images/1165979/switzerland_1721_p2_grande.jpg
https://i.pinimg.com/originals/47/be/c9/47bec90a7c79c3cfeaef6cb6c413914a.jpg
https://imgur.com/Qfwtx


texts: All LGAs for Australia in KML or KMZ format

texts: Can somebody please show me where to find a KML or KMZ file that includes the digital boundaries for all Local Government Areas in Australia?


texts: How to create a hectare raster from vector input data in QGIS

texts:
I have a vector file containing points with an attributed value (no. of employes at that location). I need to create a hectare raster with the sum of employes for each grid cell. I'm using QGIS and have so far played around with gdal_rasterize gdal_grid but with no success. I'm not even sure if I'm looking at the right tool.
Can anyone give me a hint?


texts: How to normalize the data mapped on a shapefile in GeoPandas when polygons do not have the same size

texts:
I have a shapefile containing the map of Milano that contain square polygons with different size, I mean small, medium and big squares and each square has an identifier and I want to fill these polygons with the data which is reported in another file. The Shapefile polygons are as follows:

                  id                                           geometry
0           3939_1_1  POLYGON ((504174.2673271392 5003118.268122713,...
1         3939_1_2_1  POLYGON ((508268.9938896392 5003118.268122713,...
2     3939_1_2_2_3_3  POLYGON ((511851.8796318267 5005064.825739901,...
3     3939_1_2_2_3_0  POLYGON ((511340.0388115142 5005064.825739901,...
4       3939_1_2_2_0  POLYGON ((510316.3571708892 5004415.973200838,...
5         3939_1_2_3  POLYGON ((510316.3571708892 5005713.678278963,...

This is also the data file that I have

3939_0_0    5.906552126050579   848.3214545643237
3939_0_1    5.82302081392767    1465.1476909159426
3939_0_2    11.465018361277856  742.8107082619533
3939_0_3_1  3.820090639450465   325.98139282061385
3939_0_3_3  5.1252845330079415  558.5800634528392

and I have plotted these data on Milano shapefile and I have created heatmaps but since the squares do not have the same size the data must be normalized over each polygon area. Is there any suggestion how can I normalize my data over each polygon area?

here is some part of the code that I have created these heatmaps.

import pandas as pd
import geopandas as gpd
from geopandas import GeoSeries, GeoDataFrame
mapDir = "/home/foroogh/PhD/milano-grid"
mapName = "intersection_Milano_W_GRIDIT_NEW.shp"
map_path = os.path.join(mapDir, mapName)
milano = gpd.read_file(map_path)
.
.
.
df1 = pd.read_csv(file_,sep='\t',index_col=None, header=None, names=    ['squareid','nofcalls','callduration'])
df1 = milano.set_index("id").join(df1.set_index('squareid'))
fig = plt.figure()
df1.plot(column = "callduration", cmap = "jet", figsize = (20,20))
filename3 = str(file_.split("/")[-1])+'.png'
plt.savefig(os.path.join(root_path2, os.path.basename(file_)))
plt.savefig(os.path.join(root_path2,filename3))
plt.close(' ')


texts: Converting a California Coordinate System Zone 6 value

texts:
On this page: State Plane Coordinate System - Convert, View on Google Earth I have entered:

0406 California Zone 6
X (meters) 1895486.55
Y (meters) 329169.63

Then "Calc" produces:

Degrees Minutes Seconds     30°37'18.7129", -117°20'20.1205"

Using Google Earth it lands off the west coat of Baja California rather than in San Diego County, California  near Julian as expected given the source document:


texts:
I use flex and OpenScale and I want to limit visibility just over Morocco using the coordinates, please how can I do that?!

    <os:Map id="fxmap" width="100%" height="100%" zoom="3" center="4.83212,45.75781">
    <os:Mapnik name="base"/>
    <os:DragHandler/>
    <os:ClickHandler/>
    <os:WheelHandler/>
    <os:PanZoom map="{map}" x="{fxmap.x+10}" y="{fxmap.y+10}"/>
</os:Map>

I can do that on the properties of the map ?


texts: Rasterio doesn't show the coordinates when using "imshow "

texts:
I have sentinel 2 image and I use rasterio to display it. I am using "imshow" to display it. when I display it, it doesn't show any coordinates:
src=rasterio.open('14072020.tif')

#Display the image
plt.figure(figsize=(20,10))
plt.imshow(src.read(4), cmap='gray')
plt.title('Gray Scale - Band 4')

I have ried to add it like this:
plt.imshow(src.read(4),transform=src.transform, cmap='gray')

but then I get error message:

TypeError: 'transform' must be an instance of
matplotlib.transforms.Transform, not a affine.Affine

the weird thing is that if I display 3 bands together it appears:
#RGB image
plt.figure(figsize=(20,10))
show(src.read([4,3,2]),transform=src.transform,title='Image- bands 4,3,2 IN 0-255 Values',vmin=0,vmax=30000)

Why is this happen? is there any way I can control add coordinates to the one band image?


texts: Show selected features on all zoom levels in ArcMap?

texts:
Is there a way to display selected features no matter the scale level in ArcMap 10.3? When I select by attribute and then zoom to selection I often just get a blank map since the layer is set not to display at the resulting zoom level.

To just change the layer to display no matter the zoom level is not an option, since that would result in the whole map being covered in points. The only workaround I have been able to come up with is to (1) copy the layer, (2) set it to display at all zoom levels, and (3) do the selection as a definition query. In my opinion that is to many clicks for such a simple task. Is there a better way?


texts: Is it possible to merge polygons, using ogr2ogr cmdline tools? I.e. taking a SHP file and generate a new one, where records are merged when a field has the same value?


texts: Mapnik render Error : Could not create datasource for type 'postgis'

texts:
I try to use mapnik render for a tile server. During the initialization 'renderd -f -c /usr/local/etc/renderd.conf', the following error message is displayed :

An error occurred while loading the map layer 'xxx': Could not create
datasource for type: 'postgis'  encountered during parsing of layer
'line features' in Layer at line 19 of 'XXX/test.xml'

There is no "No plugin found for type ' postgis '" displayed.

The file XXX/test.xml is the following :

<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE Map [
<!ENTITY % entities SYSTEM "inc/entities.xml.inc">
 %entities;
 ]>

<Map bgcolor="#B5D2CE" srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over">
<Style name="all_roads">
 <Rule>
  <LineSymbolizer>
   <CssParameter name="stroke">#800080</CssParameter>
   <CssParameter name="stroke-width">1.4</CssParameter>
  </LineSymbolizer>
 </Rule>
</Style>

<Layer name="line features" status="on" srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over">
 <StyleName>all_roads</StyleName>
  <Datasource>
   <Parameter name="type">postgis</Parameter>
   <Parameter name="password">XXX</Parameter>
   <Parameter name="host">localhost</Parameter>
   <Parameter name="port">5432</Parameter>
   <Parameter name="user">XXX</Parameter>
   <Parameter name="dbname">XXX</Parameter>
   <Parameter name="table">(select way from planet_osm_line where highway is not null) as roads</Parameter>
   <Parameter name="estimate_extent">false</Parameter>
   <Parameter name="extent">-20037508,-19929239,20037508,19929239</Parameter>
 </Datasource>
</Layer>
</Map>

Can you help me? Thanks


texts: MapInfo / Mapbasic: simple compile of table from six source tables

texts:
I have a project which I am reviving after a gap of some years. It involves several hundred runs within MapInfo, scripted in MapBasic.

I am stuck at a point that is so simple, it's driving me nuts.

At one point within each run, I have six tables, each of a single cell.

These tables are called:

minBNGy maxBNGx minBNGy maxBNGy medBNGx medBNGy

They are min, max and median British National Grid coordinates. Each table is mappable, with a single point.

In each case, the table's single variable has the same name as the table.

Each is a long integer, but for some (probably deliberate) reason I no longer understand, I created these as floats.

I need to compile from these six tables a single table, with the six named cells and their values on a single row. It does not need to be mappable.

I was using this:

Add Column TARGET_TABLE (NEW_VARIABLE)From SOURCE_TABLE
   Set To SOURCE_VARIABLE
   Where COL1 = COL1

... but it no longer seems to work.

In my old code, the 'map' was dropped later.

I am using MapInfo / MapBasic v.6.5 on Windows 7 Ultimate 64-bit.

I would be very grateful for any guidance.

Bruce


texts: gdalcalc.py in Python 2.7 subprocess Popen - not spawning, nor any Exceptions

texts:
from subprocess import PIPE, CalledProcessError, Popen
import subprocess
import osgeo

raster = '/home/rphillips/Desktop/tif_dir/merged.tif'
raster_dir = os.path.dirname(raster)
binary = raster_dir + 'binary.tif'
calc_list = []
calc_list = 'gdal_calc.py -A {} --outfile={} --NoDataValue=-9999 --calc "1 * (A != -9999)"'.format(raster,binary)
calc_for_subprocess = calc_list.split()
print calc_for_subprocess

try:
    df = Popen(calc_for_subprocess, stdout=PIPE)
    output, err = df.communicate()
    print output
    print err
except CalledProcessError as e:
    print (e.message)
except Exception as E:
    print (E.message)

-----------output-----------------
['gdal_calc.py', '-A', '/home/rphillips/Desktop/tif_dir/merged.tif', '--outfile=/home/rphillips/Desktop/tif_dirbinary.tif', '--NoDataValue=-9999', '--calc', '"1', '*', '(A', '!=', '-9999)"']

None

I ran the same variable
calc_for_subprocess in BASH Terminal;

rphillips@OL100314:~$ gdal_calc.py -A /home/rphillips/Desktop/tif_dir/merged.tif --outfile=/home/rphillips/Desktop/tif_dir/binary.tif --NoDataValue=-9999 --calc "1 * (A != -9999)"
/usr/bin/gdal_calc.py:266: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
myNDVs=1*numpy.logical_or(myNDVs==1, myval==myNDV[i])
0 .. 10 .. 20 .. 30 .. 40 .. 50 .. 60 .. 70 .. 80 .. 90 .. 100 - Done

The terminal command created the raster that I needed.

However I'm confused why if the numpy syntax was an issue in the Py script (shown in BASH) why wasn't the list in the Python Script even passed through the Popen constructor?
The output variable is empty.

I've had a look at the Python documentation on the subprocess module and had no luck finding anything that helped. Python 2.7.15 subprocess module documentation

On a side note: I think this reference explains the BASH gdal FutureError:-
FutureWarning: comparison to None

Update: I have also tried putting the simple string calc_list when calling the Popen constructor and had the same outcome as above.

df = popen(calc_list, stdout=PIPE)

also tried this, with no luck;

subprocess.call(calc_for_subprocess,shell=True)

but this worked but gave me a weird raster, pixel values do not correspond to numpy expression;

subprocess.call(calc_list,shell=True)


texts: Fastest way to know in which section of line is a point

texts:
I have a line and a point. I want to locate where is this point in the line according to imaginary sections I make (any way I want). The output should be in which number of section is the point (sections of the line are made in an ordered way)

I tried this as an example:

import shapely.geometry as spg
import shapely.ops as spo

line = spg.LineString([(0,0),(1,0),(1,1),(2,1),(2,2)])

# WAY TO SECTION I COME UP WITH
points_coord = []
number_sections = 4
for i in range(0,number_sections):
    point = line.interpolate(i+1)
    points_coord.extend(point.coords[:])
points = spg.MultiPoint(points_coord)

#SECTION OF LINE
geom = spo.split(line,points)

# FIND POINT IN LINE SECTIONS
p = spg.Point(0.5,0)
for index,sec_line in enumerate(geom):
    if (sec_line.intersects(new_p)):
        print("it is in section number:",index)

Is there any other way to do this faster?

Maybe I am not aware of a better shapely function that does exactly this or a more compact and faster way in basic python to loop all lines, etc.

I need to do this with a lot of points in real time and I need to be as fast as possible.


texts: Bringing information from one shapefile into another using ModelBuilder?

texts:
I'm trying to use the ModelBuilder to fasten up my work in ArcGIS but I don't know how far I can use it and what's possible. Here is my problem:

I want to analyse agriculture parcels regarding to the watersheds. I have two shapefiles: one with all my parcels and one with all my watersheds.

I now used ModelBuilder to Clip the parcels with each watershed.

This works perfectly and now I have many shapefiles with the parcels from every watershed. (for example: parcels_watershed01.shp, parcels_watershed02.shp,...)

What I would like to do, is use the information I can get from these parcels_watershed shapefiles and bring them into my origin watershed shapefile. For example: I would like to know the area of agriculture parcels inside of each watershed and bring this information into the watershed shapefile.

How can I do this in an automatic way? Is it possible to use ModelBuilder in order to extract an information from one shapefile into another? (I now opened each parcel_watershed shapefile, looked up the statistics of the area and wrote them down and had to fill them manually into the Attribute Table of the watershed shapefile... this is long and I'm sure there must be a better way, but I don’t know how.)

I'm using ArcGIS 10.1 for Desktop.


texts: GeoTIFF is solid black after QGIS Python export

texts:
I am trying to automate exporting hundreds of layers from a QGIS project in order to export them into GeoServer. I have a Python script that is exporting the TIFFs successfully but the exported tifs are always black. Any pointers?

Here is the information of the GeoTIFF using gdalinfo

Driver: GTiff/GeoTIFF
Files: m_3208548_se_16_1_20171020_20180102_a.tif
       m_3208548_se_16_1_20171020_20180102_a.tif.aux.xml
Size is 6964, 8236
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-85.062527050195669,32.312523964278661)
Pixel Size = (0.000008983152841,-0.000007594686085)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -85.0625271,  32.3125240) ( 85d 3'45.10"W, 32d18'45.09"N)
Lower Left  ( -85.0625271,  32.2499741) ( 85d 3'45.10"W, 32d14'59.91"N)
Upper Right ( -84.9999684,  32.3125240) ( 84d59'59.89"W, 32d18'45.09"N)
Lower Right ( -84.9999684,  32.2499741) ( 84d59'59.89"W, 32d14'59.91"N)
Center      ( -85.0312477,  32.2812490) ( 85d 1'52.49"W, 32d16'52.50"N)
Band 1 Block=6964x1 Type=Byte, ColorInterp=Red
  Min=0.000 Max=0.000
  Minimum=0.000, Maximum=0.000, Mean=0.000, StdDev=0.000
  Mask Flags: PER_DATASET ALPHA
  Metadata:
    STATISTICS_MAXIMUM=0
    STATISTICS_MEAN=0
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0
    STATISTICS_VALID_PERCENT=100
Band 2 Block=6964x1 Type=Byte, ColorInterp=Green
  Min=0.000 Max=0.000
  Minimum=0.000, Maximum=0.000, Mean=0.000, StdDev=0.000
  Mask Flags: PER_DATASET ALPHA
  Metadata:
    STATISTICS_MAXIMUM=0
    STATISTICS_MEAN=0
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0
    STATISTICS_VALID_PERCENT=100
Band 3 Block=6964x1 Type=Byte, ColorInterp=Blue
  Min=0.000 Max=0.000
  Minimum=0.000, Maximum=0.000, Mean=0.000, StdDev=0.000
  Mask Flags: PER_DATASET ALPHA
  Metadata:
    STATISTICS_MAXIMUM=0
    STATISTICS_MEAN=0
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0
    STATISTICS_VALID_PERCENT=100
Band 4 Block=6964x1 Type=Byte, ColorInterp=Alpha
  Minimum=0.000, Maximum=0.000, Mean=0.000, StdDev=0.000
  Metadata:
    STATISTICS_MAXIMUM=0
    STATISTICS_MEAN=0
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0
    STATISTICS_VALID_PERCENT=100

Here is what I have so far for the python script

myDir = 'C:/mydirectory'
pipe = QgsRasterPipe()
for layer in iface.mapCanvas().layers():
 print("Converting " + myDir + "/" + layer.name() + '_a.tif')
 extent = layer.extent()
 width, height = layer.width(), layer.height()
 renderer = layer.renderer()
 provider=layer.dataProvider()

 crsSrc = QgsCoordinateReferenceSystem(3857)
 crsDest = QgsCoordinateReferenceSystem(4326)

 xform = QgsCoordinateTransform(crsSrc, crsDest, QgsProject.instance())
 newExtent = xform.transformBoundingBox(extent)

 newExtent.provider.bandScale

 print(newExtent.xMaximum())

 pipe.set(provider.clone())
 pipe.set(renderer.clone())

 print(myDir + "/" + layer.name() + '_a.tif')

 file_writer = QgsRasterFileWriter(myDir + "/" + layer.name() + '_a.tif')
 file_writer.writeRaster(pipe, width, height, newExtent, crs)

I think I know what's happening. It appears that during the automated export my stats on my GeoTIFF are all being set to zero causing the output of the GeoTIFF to be blank. I updated my output using GDAL info above.

The question now is how do I make sure the stats are set during the import?


texts:
Probably there is a very simple solution to my problem - I just can't find one!
I have a point shapefile with x, y and z data - I can display it in arcscene but it is flat - the z data isn't being used! How can I direct arcscene to the z field in my attribute field so that it displays points at different levels?!


texts:
One is LandPlot and the other is Soils.

When I open the LandPlot attribute table I see I have 500 specific plots of land. When I open the Soils attribute table I see I have 5 types of soils (let's call them TypeA through TypeD). The issue is I want to select the specific Landplots (not all 500!) that contain TypeA AND TypeC soils (they can contain other soils in the plot but must include at least those two soils). I can't determine the proper process though.


texts: Single polygon shapefile to multiple polygon shapefiles

texts:
I created a single shapefile that contains multiple polygons by selection by location of the route I had and the ESRI USA counties polygon layer to give me the counties that were in the area I was interested in.

I want to know if there is a tool to use that can take that new shapefile I created of only the new polygons in my area of interest into an individual shapefile for each county.  So say there were 25 counties that crossed my route, I'd like to have a quick way to get shapefiles individually for each 25 counties.

I'm sure there is a quick script for this, but I'm not too familiar with scripts didn't know if there was quick easy way for this than exporting each individually into shapefiles.

I'm also using ArcMap 10.2


texts: Reading Raster Data with ArcGIS tools in C#

texts:
I have a raster data layer as IRasterLayer and I am using C# SDK of ArcGIS. Here is the IRasterLayer.

I want to access all the cells
I want to edit cell values.

After all i need how can i read the raster data. Where will I get the values in IRasterLayer


texts: Ranking target features by number of near features using ArcGIS Desktop?

texts:
To help decisionmaking, I want to rank housing units by the number of grocery stores nearby, limited to a 2-mile radius for each unit. I am trying to define my mapping problem and find out what ArcMap tools will help.

list the distance between target feature (house) and nearby feature (grocery) within 2 miles;
rank the target features by the number of near features within specific radii. For example, 1 store between 0 ft and 1320 ft radius, 1 store between 1320 ft and 2640 ft radius, 2 stores between 2640 ft and 5280 ft radius, and 4 stores between 5280 ft and 10560 ft radius.

Point Distance (from the Proximity toolset) lists the distances but doesn't include unique IDs or the names of my features, frustratingly. Also, it's only a table.

I am measuring the distances to multiple assets of different
FC type (points <> store locations, hospitals; and polygons <> parks).

Will I be using Euclidean Distance and Reclassify?


texts: I'm using ArcGIS Javascript API, and I need to query the GIS and show the results on the map.  Currently I have used a Query Task by specifying the layer,  the "where" statement, graphic symbol, InfoTemplate etc. and it works well.  I'm wondering however whether I should be using a FeatureLayer with a Definition Expression, rather than the Query Task.  Using a Feature Layer would enable the drawing information from the server to be used.  Also the symbol would show in the legend.  When is it better to use a Query Task rather then a Feature Layer?


texts: Calculating minimum distance between two point layers within polygons in QGIS

texts:
I have three vector layers two points and one polygon:

Centroids (POINTS FROM)

NUTS_ID | CNTR_CODE | NUTS_NAME | id

Sorting_plants_2020 (POINTS TO)

postcode | name | where | capacity | id

National_BC / countries (POLYGON)

CNTR_CODE | NUTS_NAME | BC_Balance | id

I want to calculate the minimum distance of each point to another point layer if they are located within the same country (i.e. the distance from the closest point of the closest point).
Is there any available plugin in QGIS for doing this? Or maybe some existing PyQGIS code?


texts: Publish model in toolbox as Geoprocessing Service?

texts:
I have created a model. This models takes three shapefiles and returns the result in the form of shapefile. I run this model and found it works fine and a shapefile is created and visible in the result window and I can add it to my table of contents window.

I then publish this model as a Geo-processing service and then accessed this service. The model runs successfully as the message from the server is submitted executed and successful. There are four files created in the arcgisjobs directory of the server definition, message, result and status. These are .DAT files of which the result file is largest.

However I cannot see any shapefile created neither in the result window nor in the arcgisjobs. The only thing I can see in the result window is that name of out file with message .

I added this to my table of contents and upon looking at its properties found that it has a type of inmemory feature class, however, it does not contain any data.

What do I need to do so that the output will have some data?

I am using ArcGIS 10.0.


texts: Get lat and long in two columns as floats for all points in MULTILINESTRING in PostGIS?

texts:
ST_DumpPoints(geom) gives me something like this:

("{1,1}",0101000020E61000002DEDD45CEEB362C0CE1951DA1BC24E40)
("{1,2}",0101000020E6100000C94AD826F1B362C09D4B7BC72CC24E40)...

Where I want it to look like this

seq   lat   lon

1,      61.1,  -149.1

2,      61.2,  -149.2
...

I've played around with ST_AsText and ST_AsGeoJson which gets me closer but I would like to parse it with SQL instead of the text.


texts: Add qlr group layer definition to existing QGIS map document

texts: I am able to drag and drop a .qlr group layer definition into new QGIS map projects prior to a first save, but when I drag and drop a .qlr file into a map project in QGIS after the initial save, none of the data or styles are added to the map. The names of the layers in the group appear in the browser panel, but the layers and styles are gone. Does anyone know of a solution to this problem? I am using QGIS 2.18.3.


texts:
I am trying to make these two methods variable, so the user could choose which method he wants to use. If the one with 3 classes or with 5.

    arr = arcpy.da.FeatureClassToNumPyArray(in_features, field)

    arcpy.AddField_management(in_features, "Class", "SHORT", 5, "", "","", "NULLABLE")

    methods = Method

    #METHOD THREE
    arrp = arcpy.da.FeatureClassToNumPyArray(in_features, field)
    arr = np.array(arrp,np.float)

    ##3 ranks
    p1 = np.percentile(arr, 33)  # rank = 1
    p2 = np.percentile(arr, 67)  # rank = 2
    p3 = np.percentile(arr, 100)  # rank = 3

    #use cursor to update the new rank field
    if Method[0] == "3":
        with arcpy.da.UpdateCursor(in_features , [field,'Class']) as cursor:
            for row in cursor:
                if row[0] < p1:
                    row[1] = 1  #rank 0
                elif p1 <= row[0] and row[0] < p2:
                     row[1] = 2
                else:
                     row[1] = 3
                cursor.updateRow(row)
    print Method

    #METHOD FIVE
    #arrp = arcpy.da.FeatureClassToNumPyArray(in_features, field)
    #arr = np.array(arrp,np.float)

    ##to create 5 ranks
    p1 = np.percentile(arr, 20)  # rank = 1
    p2 = np.percentile(arr, 40)  # rank = 2
    p3 = np.percentile(arr, 60)  # rank = 3
    p4 = np.percentile(arr, 80)  # rank = 4
    p5 = np.percentile(arr, 100) # rank = 5

    #use cursor to update the new rank field
    if Method[0] == "5":
        with arcpy.da.UpdateCursor(in_features , [field,'Class']) as cursor:
            for row in cursor:
                if row[0] < p1:
                    row[1] = 1  #rank 0
                elif p1 <= row[0] and row[0] < p2:
                     row[1] = 2
                elif p2 <= row[0] and row[0] < p3:
                     row[1] = 3
                elif p3 <= row[0] and row[0] < p4:
                     row[1] = 4
                else:
                     row[1] = 5
                cursor.updateRow(row)
    print Method

and adjust this:

     #parameter
    methods = [3,5]

    Method = arcpy.Parameter(
        displayName="Choose Only One Number Of Classes",
        name="method",
        datatype="Long",
        parameterType="Required",
        direction="Input",
        multiValue=False,
        enabled=True)
    Method.filter.type = "Valuelist"
    Method.filter.list = methods

The script runs but creates just classes according to Method 5.


texts: Setting parameters in ArcGIS desktop script

texts:
How would I create the parameters pixel_type and number_of bands (as found in  MosaicToNewRaster) as user inputs at the start of my own script?

ERROR 000714: Error in script BatchMosaic.
Error in executing: cmd.exe /C C:\Erosion\Scripts\BATCHM~1.PY  "C:\finalproject\Mosaic" "C:\finalproject\Mosiac_out" "mosaicme" "32_BIT_SIGNED" "1"
Failed to execute (BatchMosaic).

import arcpy, os, sys

#User parameters.
in_workspace = sys.argv[1]
out_workspace = sys.argv[2]
output_name = sys.argv[3]

def getParameterInfo(self):
    param4 = arcpy.Parameter(
        pixel_type="Input value")

    param4.filter.type = "ValueList"
    param4.filter.list = [""]

def getParameterInfo(self):
    param5 = arcpy.Parameter(
        number_of_bands="Input value")

    param5.filter.type = "ValueList"
    param5.filter.list = []

#Attempts to mosaic all raster files in root and sub directories.
rasters = []
for dirpath, dirnames, filenames in arcpy.da.Walk(in_workspace, topdown = True, datatype="RasterDataset", type = "TIF"):
    for filename in filenames:
        rasters.append(os.path.join(dirpath, filename))

#Pixel Type and Number of Bands must be manualy inserted to script.
arcpy.MosaicToNewRaster_management(';'.join(rasters), out_workspace, output_name, pixel_type = "32_BIT_UNSIGNED", number_of_bands = 1)


texts: Differences between points and lines as poles in an equal area projection

texts:
I am planning on using an equal area projection in a current project. I used this website to determine which projection I should be using.

I hadn't previously considered it, and was wondering if there are any hard set guidelines/rules about the advantages or disadvantages of poles being represented as points or lines (I will eventually be calculating total values of rasters throughout the globe).


texts: Digital Elevation Model and Shapefile not aligned in openWind

texts:
I have a digital elevation model and a shapefile containing point features for the location of wind turbines. These two layers appear to be aligned in QGIS but when I open them in openWind, they are not aligned.

I am guessing this has something to do with CRS but I am unsure how to align these layers?


texts: Convert WGS84 coordinate to Google Map coordinate system (EPSG::3785) with GeoTools and Java

texts:
Could you please point me out where can I make such a transformation:

Convert WGS84 (EPSG::4326) coordinate to Google Map coordinate system (EPSG::3785), and vice versa.


texts: Creating ArcMap custom layout that is blank?

texts:
In my work i frequently want to print draft maps on A4 or A3, and i want to make a layout template that i can apply to whatever i am doing at the time.
the problem is that the layout templates available either contain too many features (text and legend boxes etc) or are just ridiculous (as in the case of the ISO A4 or A3 templates that have data frame extents that leave a huge blank border around the map).

If i follow any of the help out there i can make a .mxd layout but this includes the data extent and map layers. i dont want to specify this as in a week my data layers may have changed or i am working in another area of the same region.... or somewhere else entirely!

How can i have a standard layout that i can apply to whatever i am working on, that gives me a page size, data frame in relation to the page and simple features like scale bar, north arrow etc?

It is getting iritating changing all this everytime i change page size.


texts:
I am using ArcGIS desktop. I want to run a tool to find all points that are close to each other within a specific distance.

What is the fastest and best tool to use?
I am using the standard arcgis with the Spatial Analyst extension.


texts: Performing viewshed in QGIS with GPS satellites?

texts:
I want to do viewshed analysis in QGIS.

My observer points are gps satellites the actual elevations are scaled down from 35000 kms to 1500m but the positions are accurate i,e elevation and azimuth.

This I am able to do in ArcGIS Desktop with not much difficulty by defining input raster, observation points(gps satellites), surfface offset (offsetB), observer elevation i,e 1500m and offsetA (1).

How do I define the similar parameters in QGIS and what shall be the search radius (default shown as 5000)?

No output is coming in QGIS.


texts: User friendly list of processing algorithms in the QGIS 3 python console

texts:
The QGIS user guide section Using processing algorithms from the console has not yet been updated for QGIS 3 and while a previously answered question will do the trick with some extra work, it is not as easy and user friendly as the old alglist()

In the user guide the following commands are presented, resulting in a user-friendly list:

import processing
processing.alglist()

This still works in the QGIS 2.18.18 python console. Using an algorithm's id one can then use e.g.

processing.alghelp("qgis:dissolve")

to learn about the algorithm's required parameters etc.

Are there any equivalent functions in QGIS 3?

Using something like print([str(e.provider().name()) + ":" + str(e.name()) for e in QgsApplication.processingRegistry().algorithms()]) gets me closer, but still not equivalent in readability to alglist().


texts: Batch uploading images into Google Earth

texts:
I am looking to visualize images on Google Earth. I have a couple thousand images that I would like to batch upload in order to verify results from a field study. I have a KMZ file with placemarks where the ID field is a image name that was taken in that location.
What I would like to be able to do is have some way of batch loading each image into the placemark (similar to the manual way of doing it by opening placemark properties and adding a local image to the placemark). I also have the images georeferenced (which I did with exiftools) if there is some easier way to upload them directly instead of using the placemark route.


texts:
Is it possible to change the projection used in OSM2PO? I used the following command to generate a SQL file which I loaded into my database:

java -jar osm2po-core-4.2.30-signed.jar prefix=at "filename.osm.pbf.pbf"

This produce the file at_2po_4pgr.sql. I then loaded this file into my postgres database:

psql -U postgres  -d routing -a -f at_2po_4pgr.sql

As this table is based on SRID:4326, I wondered if there was any way of changing it so that I can use a different projection system. The columns in the table are as follows:

  id integer NOT NULL,
  osm_id bigint,
  osm_name character varying,
  osm_source_id bigint,
  osm_target_id bigint,
  clazz integer,
  flags integer,
  source integer,
  target integer,
  km double precision,
  kmh integer,
  cost double precision,
  reverse_cost double precision,
  x1 double precision,
  y1 double precision,
  x2 double precision,
  y2 double precision,
  geom_way geometry(LineString,4326),
  CONSTRAINT pkey_at_2po_4pgr PRIMARY KEY (id)

I know that I can use ST_Transform to convertgeom_way` into a different projection but I was unsure how to update any of the other values as I presumed some were based on the initial projection.


texts:
I am trying to add and calculate a field for earthquakes selected within blast zones using the guidelines in this answer Is it possible to programmatically add calculated fields?. Can anyone tell me what I am doing wrong? Here is my code:

import processing

processing.runalg("qgis:selectbylocation",
"/Users/pgcseismolab/Desktop/qgis/moment_all.shp",
"/Users/pgcseismolab/Desktop/qgis/all_buffers.shp", ['within'], 0,0)

from PyQt4.QtCore import QVariant
from qgis.core import QgsField, QgsExpression, QgsFeature

moment_all = iface.activeLayer()

moment_all.startEditing()

blast = QgsField(‘myattr’, QVariant.Int )
moment_all.addAttribute(blast)
idx = moment_all.fieldNameIndex(‘myattr')
moment_all.commitChanges()

processing.runalg(“qgis:fieldcalculator”,
“/Users/pgcseismolab/Desktop/qgis/moment_all.shp”, “blast”, 0, 10, 3,
False, ““blast” = ‘b’”, None)

I used

execfile(r'/Users/pgcseismolab/Desktop/qgis/blast.py')

and my output error was:

    Python Console
Use iface to access QGIS API interface or Type help(iface) for more info
execfile(r'/Users/pgcseismolab/Desktop/qgis/blast.py')
Traceback (most recent call last):
File "<input>", line 1, in <module>
File "/Users/pgcseismolab/Desktop/qgis/blast.py", line 19
blast = QgsField( âmyattrâ, QVariant.Int )
SyntaxError: invalid syntax


texts:
I am trying to use Leaflet (https://leafletjs.com/) javascript library in a project built on ASP .NET MVC5.

For the sake of simplicity, let's just say I am trying to run the following code (provided by this tutorial: https://leafletjs.com/examples/crs-simple/crs-simple.html) on a view (Index.cshtml), which uses a predefined layout (_Layout.cshtml):

Index.cshtml

In the layout view, I placed the following code in the head of the HTML code, as instructed by these guidelines: https://leafletjs.com/download.html

<link rel="stylesheet" href="https://unpkg.com/leaflet@1.3.4/dist/leaflet.css" />

But this doesn't work and I can see that the library is not recognized even by Intellisense.

I have also tried downloading the package and as indicated, I added the files to the project, but having in mind the MVC project structure, I added the leaflet.css file to the Content folder and the rest of the js files to the Scripts folder. And in the layout view, I added:

<link rel="stylesheet" href="~/Content/leaflet.css" />

Now, Intellisense recognizes the library but still, when running the project the map is not shown in the view.

Am I missing something obvious here, maybe an import?

I am a beginner at using .NET framework and the ASP .NET MVC structure and I am still trying to figure out my way across this framework.

Realized I was loading from the Index view so leaflet.js was not being loaded first.

I can now see the leaflet map section but the PNG image is not being loaded. I believe leaflet is not being able to find so where should I place the image in the project structure? I tried placing it in the Content and Scripts folder but no luck.


texts: How to select GeoJSON vector feature in Leaflet?

texts:
I've been using OpenLayers for years,instead of comparing OL & Leaflet libraries,I decided to use Leaflet for the job for now.

I have a GeoJSON data layer where I show point objects only.Just as I have used in OpenLayers,when user hovers(which is supplied with leaflet) & toggle click where I should change style of vector(s) (mostly color).I cannot find this event system in leaflet.

Is there any way to do this ?


texts: How can i save and load items in canvas format in leaflet?

texts:
I've created a leaflet map and added draw control. I can draw lines, polygons, and markers but I need to add a canvas and draw my favorite drawing items like this :

I want to save these items and load them again on the map. How can I handle this?


texts: Missing connections between line and points

texts:
I have two shapefiles

line layer (Blue color)
points layer which represents houses (Green color))

I am trying to connect line with points (houses). For that, I used closestSegmentWithContext as shown in the below script of `PyQGIS.
Problem I am facing: The PyQGIS script connects only few points. There is no connecting line between points and line at some points as shown with red crosses in below image.

PyQGIS script:
roadlyr = QgsProject.instance().mapLayersByName('Line')[0]
addrlyr = QgsProject.instance().mapLayersByName('Point')[0]

#List all features
adresses = [f for f in addrlyr.getFeatures()]
roads = [f for f in roadlyr.getFeatures()]

##Create vector layer
#vl = QgsVectorLayer("LineString?crs={}" + "&field=id:integer""&field=distance:double(20,2)&index=yes".format(roadlyr.crs().authid()), "myLayer", "memory")

uri = "LineString?crs={}:"  + "&field=id:integer""&field=distance&index=yes" # Generate Conducting_equipment_mrid

vl = QgsVectorLayer(uri,
                      'Connections',
                      'memory')

provider = vl.dataProvider()

feats = []
for adress in adresses:
    a = min([road.geometry().closestSegmentWithContext(adress.geometry().asPoint()) for road in roads], key=lambda x: x[0]) #For each road find closest address
    gLine = QgsGeometry.fromPolyline([QgsPoint(adress.geometry().asPoint()), QgsPoint(a[1])])

    f = QgsFeature()
    f.setGeometry(gLine)
    f.setAttributes([adresses.index(adress),f.geometry().length()])

    feats.append(f)
    provider.addFeature(f)

QgsProject.instance().addMapLayer(vl)


texts:
I would like to calculate travel time by bicycle in postgis. Following this example from the book QGIS becoming a GIS power user, the calculation is in walking distance. How can I calculate if bicycle travel means are used.

ALTER TABLE newark_osm ADD COLUMN traveltime_min float8;
UPDATE newark_osm SET traveltime_min = length_m / 6000.0 *
60;
SELECT *
FROM pgr_drivingdistance('SELECT id, source, target,
traveltime_min as cost FROM newark_osm'::text, 1,
100000::double precision, false, false) di (seq, id1,
id2, cost)
JOIN newark_osm rd ON di.id2 = rd.id;


texts: Installing ArcGIS Desktop , QGIS and independent Python packages on same machine?

texts:
I would like to have on the same machine all my GIS software installations. Both ArcGIS Desktop and QGIS come with their own Python installations but I want another independent Python where I can play around and update its packages without messing up ArcGIS and QGIS Pythons. I understand that just installing the independent Python into another folder should suffice.

However, as I use GDAL installed from GISinternals files and pip to update packages, what should I be careful with in order to install GDAL Python bindings and pointing to the correct Python version?

I think I recall, when installing Python bindings it asks for the instance of Python you are to install it for, is that correct?

Also, when using pip for installing packages, how can I make sure I am using the independent Python's pip instead of, for example, ArcGIS's?


texts:
Using ArcMap 10.1 I Merged multiple .shp layers (each with multiple polylines) into one .shp

The resulting layer will not draw, I get the following error message: One or more layers failed to draw: [LAYER NAME] There was a problem drawing the shapefile. Drawing aborted.

I have re-done the merge but have the same problem.  This is a very large layer - is there a size limit? Could it be taking too much memory to draw?

I don't need to see the layer but want to run intersect with it (so that I can calculate the length of the lines in multiple buffers) and that is returning a blank output - I presume because it is failing to draw.


texts: Is it possible to make a label with underline and arrow from label to the point in GeoServer's SLD?

texts: It is commonly a problem when you need all points to be labelled, but due to high density of points it is difficult to understand the ties between points and labels. Arrows would be a way, but I haven't found anything similar even in SLD extensions.


texts: How to recover deleted records from a geodatabase?

texts: A coworker of mine made a copy of a personal geodatabase feature class, and deleted a subset of records from the attribute table.  He was not aware that the original database was still linked to his copy, and therefore he deleted records from the original database.  Is there any way to recover these records?


texts: In WMS service Geographic Coordinate System changed from NAD 1983 to WGS 1984

texts: In our project MXD it's Geographic Coordinate System is "NAD1983" and Projected Coordinate system is "NAD_1983_BLM_Zone_11N", and I published Map and WMS service from this MXD, while publishing in capabilities I checked mark on both Mapping  and WMS. After publishing services I opened both services in separate blank MXDs for Map service Geographic Coordinate System and Projected Coordinate system are same as source MXD but in WMS service Geographic Coordinate System changed to WGS_1984 and for Projection it showing supported coordinate system are WGS_1984, NAD_1983_BLM_Zone_11N. I am not getting why in WMS service Geographic Coordinate System changed to WGS_1984? Although in source MXD layer and dataframe has same GCS and PCS.


texts: How could one use a GIS (or other appropriate system) to calculate the change in mean ocean levels across the world if all the land below sea level were flooded (Death Valley, Caspian Depression, Afar Depression for example)?


texts: How to split/divide linestrings of larger lengths in postgis

texts:
I want to split/divide linestrings in my data. I have roads data and I want to have every road segment having length not more than 250 meter.

Check the shape length of this particular line feature. I want every road length less than 250 meter. Obviously this particular road must be divided in several parts.

If ST_Line_Substring  is the solution, kindly elaborate. I can't figure out how to apply this in my case


texts: Geoserver - Geofence How to change catalog mode?

texts:
According to the docs it seems to be possible to change the catalog mode for a layer using this Geofence plugin.

But I can't find an option, button or something like that, that could do this in GeoServer page.

How it is possible to change the catalog mode for a layer in Geofence? Is it only possible via REST?

Changing the catalog mode in Layer Security section seems to have no effect at all.

I'm using GeoServer 2.11-RC1


texts: ArcGIS Geometric Network - Identify Indeterminate Flow

texts: We have a geometric network with sewer pipes that have their flow direction based on digitized direction.  Some of the pipes, however, have indeterminate flow.  We would like to find all the pipes that have indeterminate flow.  We can do this by setting the flow direction arrows and identifying these pipes visually/manually, however, since this is a rather large network, is there a way/tool that can identify all the pipes that have indeterminate flow in a simpler, en masse method?  We are using ArcGIS 10.2.2.


texts: GeoTools (Java) read properties from ESRI shapefile at specific coordinates

texts:
This is my first time using GeoTools and I have a 3 simple things I need to get from .shp files which I'm reading with my program.

I need to know in specific coordinates (latitude/longitude).

Example: lat=58.85 lon=5.26 for 3 specific things:

Exclusive Economic Zone (EEZ)
FAO area
ICES statistical rectangle

I know that all of these properties are there, because tested the whole shapefile and it gave my these results:

EEZ<EEZ id=eez.232>=German Exclusive Economic Zone
EEZ<EEZ id=eez.233>=Azerbaijanis Exclusive Economic Zone
EEZ<EEZ id=eez.234>=Kazakh Exclusive Economic Zone

And same is for ICES rectangles:

ICESNAME<ICESNAME id=ices.11002>=26M8
ICESNAME<ICESNAME id=ices.11003>=27M8

It reads the whole shapefile, but I don't know how to read these 3 properties from specific coordinate.

My program code:

File file = new File(shapeDirectory + "ices.shp");

Map<String, Object> map = new HashMap<>();
map.put("url", file.toURI().toURL());

DataStore dataStore = DataStoreFinder.getDataStore(map);
String typeName = dataStore.getTypeNames()[0];

FeatureSource<SimpleFeatureType, SimpleFeature> source = dataStore.getFeatureSource(typeName);
Filter filter = Filter.INCLUDE; // ECQL.toFilter("BBOX(THE_GEOM, 10,20,30,40)")

SimpleFeatureType schema = source.getSchema();
String geomType = schema.getGeometryDescriptor().getType().getBinding().getName();
System.out.println(geomType);

    FeatureCollection<SimpleFeatureType, SimpleFeature> collection = source.getFeatures(filter);
    try (FeatureIterator<SimpleFeature> features = collection.features()) {
        while (features.hasNext()) {
            SimpleFeature feature = features.next();
            System.out.print(feature.getID());
            Property property = feature.getProperty("ICESNAME");
            System.out.println(property);
        }
    }
}

How to read all these 3 attributes (from 3 different shape files (fao/ices/eez)) at specific coordinates?
The coordinate input will be decimal, but if someone knows, please provide me the degrees/minutes/seconds example too.


texts: Reading L1B ASTER in GRASS using r.in.aster

texts:
I'm having a problem getting an ASTER L1B to be read in GRASS using r.in.aster.
The output the command shows this:

  (Wed Apr  8 17:49:03 2015) Command finished (1 sec)
  (Wed Apr  8 17:49:37 2015)
  r.in.aster
  input=/home/gobel/Documents/Images/Aster/AST_L1B_00301082001131449_20110622184108_17885.hdf
  proctype=L1B band=1 output=test Georeferencing aster image ...
  Creating output file that is 6133P x 5675L. Processing input file
  HDF4_EOS:EOS_SWATH:/home/gobel/Documents/Images/Aster/AST_L1B_00301082001131449_20110622184108_17885.hdf:VNIR_Swath:ImageData1.
  0...10...20 ...30...40. ..50...60. ..70...80.. .90...100 - done. Importing into GRASS ... Traceback (most recent call last):   File
  "/usr/lib/grass70/scripts/r.in.aster", line 168, in
      main()   File "/usr/lib/grass70/scripts/r.in.aster", line 129, in main
      import_aster(proj, srcfile, tempfile, band)   File "/usr/lib/grass70/scripts/r.in.aster", line 160, in import_aster
      outfile = "%s.%s" % (output, band) NameError: global name 'output' is not defined (Wed Apr  8 17:49:38 2015) Command finished (1 sec)


texts:
I am labeling my polygons with text-placement: interior; which places the label at the center of each polygon.
When I zoom in so that the center of the polygon with the label is not inside the map extent anymore I still want the polygon to be labeled. Ideally again in the center of the intersection between the map extent and the original polygon geometry.
Any idea how to accomplish that?


texts:
I am getting an extension error when I try and intersect some points and polygons.

Here is a photo of what I am trying to intersect:

Any help is appreciated.


texts: ST_MakeLine not returning proper latitude/longitude values


texts:
I am on Ubuntu 13.04. I get the following error message -

    Traceback (most recent call last):
      File "analyse.py", line 1, in <module>
        from log import shelve
      File "/home/shubham/SMART/TaxiData/log.py", line 27, in <module>
        from demo import *
      File "/home/shubham/SMART/zones/demo.py", line 5, in <module>
        from qgis.core import *
    ImportError: No module named qgis.core

Actually, everything was working fine till today morning. I guess this might be due to a package update.

The qgis version I am using is 1.8.0-2~raring1. The remote repository is http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu/ raring/main amd64 Packages.

I tried looking around on Google but my search was fruitless. So, I will really appreciate any help or pointers you guys can give :)

Thanks.


texts:
I want to write a binary raster to disk. I have the numpy array whose dtype is:

dtype('bool')

I try to open a GTiff for writing:

img_output = rasterio.open(
            "../../binary_output.tif",
            'w',
            driver='GTiff',
            nodata=nd,
            height=self.raster.height,
            width = self.raster.width,
            count=1,
            dtype = binary_raster.dtype,
            crs=self.raster.crs,
            transform=self.raster.transform,)

But I end up with:

TypeError: invalid dtype: 'bool'

I have also tried things like:

dtype = bool
dtype = 'bool'
dtype = "bool"

What is the problem with this? According to the docs, dtype can be any numpy dtype. How do I write a binary raster?


texts:
I am using the ST_LineLocatePoint function from PostgreSQL. According to what I have read, the description for this function is as follows:

  Returns a float between 0 and 1 representing the location of the closest point on LineString to the given Point, as a fraction of total 2D line length.

I am using a 3D LineString and I want to get the 3D distance to a point.

Is there an equivalent 3D function for ST_LineLocatePoint?


texts: is it possible to store empty geometries in an Oracle Database?

texts:
I'm trying to fix a GeoTools issue but I'm running into problems putting empty geometries into my Oracle database. Basically executing:

select SDO_UTIL.FROM_WKTGEOMETRY('POINT EMPTY') from PARISH_BOUNDARIES_BUP;

gives me the ever so helpful error message:

  ORA-29532: Java call terminated by uncaught Java exception:
  java.lang.RuntimeException ORA-06512: at "MDSYS.SDO_UTIL", line 137
  29532. 00000 -  "Java call terminated by uncaught Java exception: %s"
  *Cause:    A Java exception or error was signaled and could not be
             resolved by the Java code.
  *Action:   Modify Java code, if this behavior is not intended.

if I test in eclipse with the actual code I get:

  java.io.IOException: Error inserting features     at
  org.geotools.jdbc.JDBCDataStore.insert(JDBCDataStore.java:1697)   at
  org.geotools.jdbc.JDBCInsertFeatureWriter.flush(JDBCInsertFeatureWriter.java:125)
    at
  org.geotools.jdbc.JDBCInsertFeatureWriter.close(JDBCInsertFeatureWriter.java:152)
    at
  org.geotools.data.InProcessLockingManager$1.close(InProcessLockingManager.java:350)
    at
  org.geotools.data.store.ContentFeatureStore.addFeatures(ContentFeatureStore.java:263)
    at
  org.geotools.data.oracle.OracleGeometryOnlineTest.testInsertEmptyGeometry(OracleGeometryOnlineTest.java:96)
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)  at
  sun.reflect.NativeMethodAccessorImpl.invoke(Unknown Source)   at
  sun.reflect.DelegatingMethodAccessorImpl.invoke(Unknown Source)   at
  java.lang.reflect.Method.invoke(Unknown Source)   at
  junit.framework.TestCase.runTest(TestCase.java:176)   at
  junit.framework.TestCase.runBare(TestCase.java:141)   at
  junit.framework.TestResult$1.protect(TestResult.java:122)     at
  junit.framework.TestResult.runProtected(TestResult.java:142)  at
  junit.framework.TestResult.run(TestResult.java:125)   at
  junit.framework.TestCase.run(TestCase.java:129)   at
  org.geotools.test.OnlineTestCase.run(OnlineTestCase.java:123)     at
  junit.framework.TestSuite.runTest(TestSuite.java:255)     at
  junit.framework.TestSuite.run(TestSuite.java:250)     at
  org.junit.internal.runners.JUnit38ClassRunner.run(JUnit38ClassRunner.java:84)
    at
  org.eclipse.jdt.internal.junit4.runner.JUnit4TestReference.run(JUnit4TestReference.java:86)
    at
  org.eclipse.jdt.internal.junit.runner.TestExecution.run(TestExecution.java:38)
    at
  org.eclipse.jdt.internal.junit.runner.RemoteTestRunner.runTests(RemoteTestRunner.java:459)
    at
  org.eclipse.jdt.internal.junit.runner.RemoteTestRunner.runTests(RemoteTestRunner.java:675)
    at
  org.eclipse.jdt.internal.junit.runner.RemoteTestRunner.run(RemoteTestRunner.java:382)
    at
  org.eclipse.jdt.internal.junit.runner.RemoteTestRunner.main(RemoteTestRunner.java:192)
  Caused by: java.sql.BatchUpdateException: ORA-29875: failed in the
  execution of the ODCIINDEXINSERT routine ORA-13354: incorrect offset
  in ELEM_INFO_ARRAY ORA-06512: at "MDSYS.SDO_INDEX_METHOD_10I", line
  623 ORA-06512: at "MDSYS.SDO_INDEX_METHOD_10I", line 227

  at
  oracle.jdbc.driver.OraclePreparedStatement.executeBatch(OraclePreparedStatement.java:12296)
    at
  oracle.jdbc.driver.OracleStatementWrapper.executeBatch(OracleStatementWrapper.java:246)
    at
  org.apache.commons.dbcp.DelegatingStatement.executeBatch(DelegatingStatement.java:297)
    at
  org.apache.commons.dbcp.DelegatingStatement.executeBatch(DelegatingStatement.java:297)
    at
  org.apache.commons.dbcp.DelegatingStatement.executeBatch(DelegatingStatement.java:297)
    at org.geotools.jdbc.JDBCDataStore.insertPS(JDBCDataStore.java:1762)
    at org.geotools.jdbc.JDBCDataStore.insert(JDBCDataStore.java:1683)
    ... 25 more


texts:
I ran the tool Spatial Join(Analysis),
target and joined features are line feature class,
the join operation is Join one to one,
I set a merge rule(sum) in one of my fields(the type of the field is double)
and the match option is share a line segment with

I get ''error 000426 out of memory''.
Is there any way to avoid it?
I ran the tool with only a few lines and it works fine.


texts:
I was wondering if there is any way to protect fields so that users cannot delete or edit them. With that being said, a 'superuser' would need to be able to edit and delete at will.

working in python with arcpy; arcgis 10.3


texts: How to display indoor map using OpenLayers?

texts:
I created my osm file in JOSM, and i imported it to my database by osm2pgsql tool.
I created new layer in my osm.xml file, which selecting data with query:

select way from planet_osm_line,
(select parts from planet_osm_rels WHERE (((hstore(tags)->'level'))  = '0')) p
where osm_id = any (p.parts))

when i run my python script to generate png file i got an error:

  RuntimeError: Postgis Plugin: PostGIS: unable to query the layer
  extent of table 'planet_osm_rels' because we cannot determine the
  geometry field name. Please provide either an 'extent' parameter to
  skip this query, a 'geometry_field' and/or 'geometry_table' parameter,
  or add a record to the 'geometry_columns' for your table.

is any one can tell me what i'm doing wrong and how i can fix that?


texts:
I like to use old version of FME where only Inspector are used because it is more steady. Visual Preview is too hasty and useless for me because it will always change map view when selecting different FME -transformer.

How to disable Visual Preview in FME Workbench 2019?

There is document from Safe about Visual Preview.

http://docs.safe.com/fme/2019.0/html/FME_Desktop_Documentation/FME_Workbench/Workbench/Visual-Preview.htm


texts: Publishing dynamic raster data repository using ArcGIS server?

texts:
I am trying to publish raster datasets using ArcGIS server. The dataset grows every 10 days or so and I want the new data to appear in the same service and be able to select it using simple layer selection. Is there any way I can achieve this without destroying the existing service and re-creating it?

I want to know if there is a way I can publish the folder or geo-database containing the raster datasets once in a way where the changes in that folder/geo-database is reflected in the service.


texts:
I am creating a Java web service that receives GeoJSON.

Does anyone know a java library to manage GeoJSON?


texts:
I've just installed QGIS 2.0 on my Windows 7 computer. I'm trying to open a raster and the system reports this message:

Cannot open TIFF file due to missing codec.

How can I solve it?


texts: Creating a multiversioned view with ArcSDE and SQL Server 2008 R2

texts:
I am using the administration command sdetable -o create_mv_view to create a multiversioned view on a table in SQL Server 2008 R2.

Everything appears to be working correctly with one exception:  The SDE_get_view_state() scalar function that was included with ArcSDE will only return data for version with name set to 'DEFAULT' and owner set to 'sde'.   This is hard-coded in the function.

I have some clients who are using the DBO schema for ArcSDE instead of the SDE schema, and therefore they will not have any version data that meets the criteria.

Does anyone know why the function is checking for the schema owner to be set to 'sde' and more importantly, what is the best way to get this to work for the DEFAULT layer regardless of the owner?


texts: How to build a model creating shp polygonal objects using for that the parameters of xls file?

texts:
Could somebody advise which instrument I can use in  ModelBuilder (ArcGIS10) for build a model which will create shp polygonal objects from a table (created previously in Excel) with necessary parameters (that include coordinates of latitude and longitude and distances).

Some details.
I need to create many polygonal objects by 4 steps:

Create 2 circles with the same center (with coordinates X1, Y2 – they will be in the table) and radiuses Rmax and Rmin (they will be in the table);
Create a new polygonal object (“bagel”) = Circle (max) – Circle(min);
Create a polygon triangle by coordinates of 3 their vertexes (all coordinates will be in the table);
Find a resulting polygonal object (as a shapefile) that will be equal a space of intersection of “bagel” and that triangle.


texts: Python pyproj converting from wgs84 to Israel Grid IG 05/12 EPSG:6991

texts:
I am trying to convert from lat/long to xy in Israel Grid 05/12:

links https://epsg.io/6991 (https://epsg.io/6991.gml)

I tried this :
p = pyproj.Proj("+init=ESPG:6991") using python package pyproj but it is not defined yet.

File "_proj.pyx", line 84, in _proj.Proj.__cinit__ (_proj.c:1170)
RuntimeError: no system list, errno: 2

What is the proj code or definition string to use instead (perhaps in the format like the previous network?

+proj=tmerc +lat_0=31.73439361111111 +lon_0=35.20451694444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +towgs84=-48,55,52,0,0,0,0 +units=m +no_defs ?


texts: Failed to Generate Incremental ID with Following Python Code in the Code Block. What's the mistake in my code?

texts:
I am trying to add a field with incremental ID for a shapefile with certain feature classes. I know this can be done by using Sort tool and then use the Add Incremental ID field tool. But since I'm learning Python, and this seems like not too difficult a problem to tackle with Python, I want to use Python to handle this to enhance my Python skills.

Below is the code I put in the Code Block:

def letsdothis():
    rec = 0
    vStart = 1
    vDelta = 1
    if rec == 0:
        rec += vStart
    else:
        rec += vDelta
    return rec

And the function I use of course is letsdothis()

To show it better:

My reasoning is this:

The first one will be rec += 1, then rec will equal 1. Then the next one, since rec != 0. it will be (the new rec value) += 1, which will be 2. And so on and so on.

But the results are all 1. My understanding is that for every row cell/feature class, the function will be called from scratch therefore rec will equal 0 instead hold the value that's calculated in the previous feature class/row.

Is this the problem with my code? I am very new to writing my own Python code in GIS. I've only used other's code and tweaked them to fit my purpose before so I don't know how to really do this.

If I can make rec == 1 in the first row, how shall I make that value pass on to the nest row?


texts: Setting start point of line drawn with SketchViewModel in ArcGIS API for JavaScript

texts:
I'd like to be able to programmatically set the start point of a new line drawn using SketchViewModel.create('polyline') so that the new line starts from the endpoint of a previously drawn line.  I'm not aware of any way to do this with the 4.x API nor have I had any luck simulating a map click to fake it in (which is what I did w/ 3.x Draw to accomplish this).
Example of 3.x solution (but would prefer something more direct if possible):
import Draw = require('esri/toolbars/draw');
import EsriMap = require('esri/map');
import Point = require('esri/geometry/Point');
import webMercatorUtils = require('esri/geometry/webMercatorUtils');
import screenUtils = require('esri/geometry/screenUtils');
import { DrawOptions } from 'esri';

export class CustomDraw extends Draw {
  constructor(public map: EsriMap, options?: DrawOptions) {
    super(map, options);
  }

  /**
   * Programatically add another point to the drawing.
   * @param point
   */
  addPoint(point: Point) {
    // can't find another way to do this other than simulating a map click at the desired location
    // which Draw will handle and add as if user had clicked there.
    let mapPoint = <Point>webMercatorUtils.geographicToWebMercator(point);
    let screenPoint = screenUtils.toScreenPoint(this.map.extent, this.map.width, this.map.height, mapPoint);
    (<any>this.map).emit("click", { bubbles: true, cancelable: true, screenPoint: screenPoint, mapPoint: mapPoint });
  }
}