texts: Trouble writing .tif from memory directly into S3 Bucket

texts:
I've used question on the topic as initial help and structured my code after it. I'm warping a raster file using these kwargs
rast = gdal.Open(/vsicurl/raster_url)
wkt = rast.GetProjection()
outbounds = [-180.00, -90.00, 180.00, 90.00]
warp_kwargs = {
                'format': 'vrt',
                'srcSRS': wkt,
                'dstSRS': 'EPSG:4326',
                'outputBounds': outbounds,
                'options':['SOURCE_EXTRA=1000','CENTER_LONG=0']
            }

gdal.Warp('/vsimem/rast_1.tif', rast, **warp_kwargs)

I then used this previous question to try to upload the file from /vsimem/ directly into an s3 bucket. I changed a couple things, so my code is as follows.
(I use fp to mean filepath)
def upload_to_s3(mem_fp, s3_fp, bucket_name):
    f = gdal.VSIFOpenL(mem_fp, 'rb')
    gdal.VSIFSeekL(f, 0, 2)  # seek to end
    size = gdal.VSIFTellL(f)
    gdal.VSIFSeekL(f, 0, 0)  # seek to beginning
    data = gdal.VSIFReadL(1, size, f)
    gdal.VSIFCloseL(f)
    # Upload the raw data to s3
    s3 = boto3.client('s3', aws_access_key_id=ACCESS_KEY,
                      aws_secret_access_key=SECRET_KEY)
    s3.put_object(Key=s3_fp, Bucket=bucket_name, Body=data, ContentLength=size)

In this case size is the size of the raster and data is b'<VRTDataset rasterXSize="1440" rasterYSize="720" subClass="VRTWarpedDataset">. So this process does upload a file to s3, but it's only 3 kb when it should be around 8 mb. And when I try to download the file from s3 and plot it in QGIS, it errors out saying
ERROR 4: : No such file or directory
ERROR 5: GDALWarpOptions.Validate(): hSrcDS is not set.
gdalinfo failed - unable to open 'test.tif'.

I've tried using both s3.put_object() and s3.upload_file() and neither of them work. I've also tried using an s3.Object.put() which didn't work. I've tried putting the raster in /vsimem/temp which also didn't work. Not sure what else to try here.


texts:
I've an animated torque visualization in CartoDB. The Time Column of animation is date type: 2015-07-03T00:00:00Z
The visualization is ok, but in the mapview the timeslider write the steps like '07/03/2015' (month/day/year). And I want the timeslider write '3 JUL' (day monthshortname).

How I can do this?


texts: Cannot read property 'decode' of undefined in lrm-google

texts:
I am following this doc Lrm-google to use google direction api with leaflet routing machine. I haven't changed much of my code except replacing geocoder: L.Control.Geocoder.nominatim(), with router: L.Routing.google(), my code works when i switch back to geocoder: L.Control.Geocoder.nominatim() - the routes appears. But when I use L.Routing.google(), I get this error in my console - Uncaught TypeError: Cannot read property 'decode' of undefined

    this.control = L.Routing.control({
      waypoints: [
        L.latLng(myLocation.lat, myLocation.lon),
        L.latLng(closest_latlng.lat, closest_latlng.lng)
      ],
      // geocoder: L.Control.Geocoder.nominatim(),
      router: L.Routing.google(),
      reverseWaypoints: true,
      showAlternatives: true,
      lineOptions: {
        styles: [
          { color: "black", opacity: 0.15, weight: 9 },
          { color: "white", opacity: 0.8, weight: 6 },
          { color: "steelblue", opacity: 1, weight: 4 }
        ]
      },
      altLineOptions: {
        styles: [
          { color: "black", opacity: 0.15, weight: 9 },
          { color: "white", opacity: 0.8, weight: 6 },
          { color: "hotpink", opacity: 1, weight: 4 }
        ]
      },
      missingRouteStyles: {
        styles: [
          { color: "black", opacity: 0.15, weight: 7 },
          { color: "white", opacity: 0.6, weight: 4 },
          { color: "gray", opacity: 0.8, weight: 2, dashArray: "7,12" }
        ]
      }
    }).addTo(map);


texts:
I have a shapefile with 100k points. I want to create an atlas that would only show set of points having the same specific attribute.

gid;geom;id_coupe;symb
3572;"POINT(1.79008147341714 5.93)";"alz_02_2";"Céramique"
3573;"POINT(27.1729542533557 5.99)";"alz_02_3";"Grès"
3686;"POINT(27.1356011546138 5.91)";"alz_02_3";"Céramique"
2264;"POINT(7.64830426499203 6.2)";"alz_01_0";"Lithique"
2264;"POINT(2.35880324738061 6.2)";"alz_01_1";"Lithique"
3574;"POINT(28.1625165379164 5.97)";"alz_02_3";"Céramique"
3573;"POINT(0.609316406683922 5.99)";"alz_02_2";"Grès"
3575;"POINT(28.4974037106703 5.93)";"alz_02_3";"Céramique"
3686;"POINT(3.07935785287035 5.91)";"alz_02_2";"Céramique"
2299;"POINT(2.02291014945037 6.05)";"alz_01_1";"Lithique"
2302;"POINT(7.6597310840181 6.07)";"alz_01_0";"Lithique"
2299;"POINT(7.21508858373733 6.05)";"alz_01_0";"Lithique"
2302;"POINT(2.50870130550044 6.07)";"alz_01_1";"Lithique"
2300;"POINT(2.19214124238729 6.04)";"alz_01_1";"Lithique"
3572;"POINT(27.5621763726535 5.93)";"alz_02_3";"Céramique"
2301;"POINT(2.3191856638636 6.04)";"alz_01_1";"Lithique"
2300;"POINT(7.29670253603475 6.04)";"alz_01_0";"Lithique"
2301;"POINT(7.60792547014961 6.04)";"alz_01_0";"Lithique"
3575;"POINT(0.998974763657425 5.93)";"alz_02_2";"Céramique"
3574;"POINT(0 5.97)";"alz_02_2";"Céramique"
...

So in an atlas, the purpose is to display only points depending on id_coupe, so to speak one map with only Alz_01_0 points, one with alz_01_0 points, etc. but I have to keep all the lines (no aggregate). In this way I can use symbology with symb attribute.

I see things like this.

Apply a symbology using symb attribute to all points. "Céramique" -> red point, "Lithique" -> blue point,...
Find a way to automatically create maps with all points having id_coupe = 'alz_01_0' for the first one, with id_coupe = 'alz_01_1' for the second one, all points with id_coupe = 'alz_02_1' for the third one, up to the last id_coupe so to speak the 299th one.
Of course the symbology has to be kept for all maps.

Obviously, this seems to be a challenge unless being a python Jedi (?) I am not :/


texts: Appending fields from external tables to feature class using ArcPy?

texts:
I have some data in a relational database structure within a geodatabase i.e a feature class with several lookup tables. I need to flatten this data so that certain fields from the lookup tables are appended to the feature class attribute table. A join will not work, neither will using domains and a query table can only do an inner join which is no good.

I know this can be done in python but where do I begin?

I am using ArcGIS 10.0 and the dataset is large - 1.2 million records.


texts: Using lasmerge in QGIS 2.8.3 is not working: "...instance has no attribute"

texts:
Every time I try running lasmerge tool from LASTools within QGIS, I get the following error message:

  Algorithm lasmerge starting...
  lasmerge instance has no attribute 'addParametersApplyFileSourceIdsCommands' See log for more details

I have already installed the 2.8 LASTools fixed provided by the QGIS guide.
I have checked nothing above [input LAS/LAZ file] nor have I included anything in the additional command line parameters.

I am trying to merge 4 different LAS files.

Can someone help me understand what I am doing wrong?


texts: How to create animation from raster data?

texts:
I have some NDVI and precipitation/temperature images at the monthly time step from Jan, 2005 to Dec, 2008.

I need to create animation image to show the temporal change.

How can I do this?

Is it possible to use ArcGIS, ERDAS, ENVI?
If so, please tell me the processes step by step?
If not, what software and tool can be used to do this?


texts: Transforming descriptive CSV data to known locations in GIS using Python

texts:
I have to transfer CSV data that describes a FROM and TO segment into spatial data. The segment descriptions do not follow any logical convention, but rather describe road junctions, places, or road features (4-lane section). I have decided to use a Python script to iterate through the CSV data, extract the Segment ID, highway #, FROM, and TO fields into a dictionary {ID:segHwy,segFrom,segTo}. This works fine.

The dataset that I am using to then locate these segments is an intersection of 3 datasets in order to capture all possible road names for a match. I need to now figure out how to search through this large dataset and compare it against my dictionary to find matches of all three values in my dictionary.
If both the highway number and the FROM field finds a road match, then I populate a new FROM field in the spatial dataset with the segment ID and 'from'. Same with TO.

In the end, I'm hoping to have a spatial dataset with two populated fields (Seg_From and Seg_To) with the Segment ID from CSV and whether this was a 'from' or 'to' location.
I'm mostly stuck on how to extract key from value and also search through 3 different road name fields in my spatial dataset.
I just have the dictionary so far, am using ArcPy....any help?

with arcpy.da.SearchCursor(segData,[segID, segHwy, segFrom, segTo]) as segCursor:
for row in segCursor:
    if row[0] not in segDict:
        segDict[row[0]] = row[1:]

Return: {1: ('H35', 'Jct H16', 'end Trestle'), 2: ('H16', 'Nouch Rd', 'Lake Kathlyn Rd W'),....}


texts: Where can I obtain 3D roadside features for UK?

texts:
I am doing some work on sight distances and I need a map or whatever of the UK that contains roadside features (road signs, trees etc) as well as the road features (road centreline, etc) so I can convert to an appropriate format for use in ArcGIS. I will need to eventually write a programme such that based on the required sight distance I can create a line that will intersect with these roadside features and know if they obstruct visibility.

So first I need a map of the UK which shows these features such that they can be understood by ArcGIS. If anyone knows anything please let me know.


texts: Previously I loaded segp1 files into google earth using a earthpoints program, and then I loaded the google earth kmz file into QGIS project, I want to find a better way? I think there is a GDAL plugin but I am very new to GIS and don't know where to add plugin?


texts: Arcpy geoprocessing script running with default inputs

texts:
I have a pretty simple geoprocessing script. My javascript is supposed to be passing it a list of layers to query. Something like this for the Javascript:

var params = { "Layers": "foo,bar,hello,world"};
var gp = new Geoprocessor(*myURL*);
gp.submitJob(params, gpSuccess, gpStatus, gpFailed);

The Geoprocessor just explodes the list based on the comma, and reads the layers from the database. Pretty straight forward. I'm reading the inputs in as a text string, like so:

layers = GetParameterAsText(0)

My code works fine in my own testing, the problem is when I have this running as a service, I can't seem to pass it the correct values. It always runs based off the default values (which end up being whatever I supplied when I did my test run in Arc Catalog). I have "layers" defined as an input variable for the model, so I'm not sure why it's not getting my input. Checking the job status verifies that it is running off the default values, and not what I am attempting to supply as an argument.

The geoprocessing script can be accessed through the REST interface here.

The source code can be found here.


texts:
I am trying to rasterize a vector with the QGIS 3.4 rasterize tool (not saga or grass tools). I prefer using this tool cause it is much faster. The values which i would like to burn  are found in the "cost_value" field(there is only one value in this case which is 2934). The no data value should get -9999. After i run the tool all of the values from the raster get the values from "cost_value" field. Simply said: all the pixels have the value from "cost_value" and none of them get -9999 even though i am rasterizing only one line geometry.

My settings:


texts: Displaying buses from a transport network (newbie)

texts:
First of all, I'm noob in this subject so please give me some advice.

I'd like to develop something like http://traintimes.org.uk/map/london-buses/#73 but with my own city.
Unfortunately our MTA doesnt provide TFL-like feeds, but I can get some data in realtime such as:
"Line 1, Stop 1, next bus is at 500m, another at 1000m."

1- I'd like to know how do I store and display lines and stops (using leafleat or another tool you could recommend)

2- How do I show buses, if I have the distance in meters the bus from any stop.

The page will reload manually at first.
Thanks!!


texts: Update rows based on other rows as script tool?

texts:
I am trying to convert python script to tool. Tool would update row values based on other row values. For example, if in row One values are Water, Rivers, Streams, in row Two, updated values would be "Yes", "No", "Maybe".

Using this script in this post Using UpdateCursor to update a field based on other field values, I have to managed to write my code to update multiple rows.

This is my working code, for python and one row:

import arcpy

fc = r'D:\D_Temp\Folder\Data.shp'
with arcpy.da.UpdateCursor(fc, ["Description", "New"]) as cursor:
    for row in cursor:
        if row[0] == "Water":
            row[1] = "Yes"
            cursor.updateRow(row)
        if row[0] == "Rivers":
            row[1] = "No"
            cursor.updateRow(row)
        if row[0] == "Streams":
            row[1] = "Maybe"
            cursor.updateRow(row)

EDITED (with part of the code FelixIP suggested):

When I'm trying to convert the script as python tool script, I do it like this. And I also put correct parameters (feature class as input data, field for original data, and field for derived data):

import arcpy
fc = arcpy.GetParameterAsText(0)
fieldList = arcpy.GetParameterAsText(1)
fieldList = fieldList.split(";")
with arcpy.da.UpdateCursor(fc, fieldList) as cursor:
    for row in cursor:
       if row[0]  == "Water":
           row[1] = "Yes"
           cursor.updateRow(row)
       if row[0]  == "Rivers":
           row[1] = "No"
           cursor.updateRow(row)
       if row[0]  == "Streams":
           row[1] = "Maybe"
           cursor.updateRow(row)

but no matter the parameters i define, I end up with this error:

  Traceback (most recent call last):
    File "C:\Users\Name\Desktop\Script.py", line 7, in
      row1 = "Yes"
  IndexError: list assignment index out of range.
  Failed to execute (NewScript).

I think I am missing some sort of list for rows, but I am not sure.


texts: How to Access the Current Workspace Using an ArcObjects Add-in for ArcMap?

texts:
Accessing current workspace using an ArcObjects add-in for ArcMap (Table of Contents has layer groups).

I am trying to retrieve the current WorkSpace using ArcObjects for an ArcMap add-in in C#.
Any help would be appreciated. I've tried this:

IWorkspace ws = ((IDataset)ArcMap.Document.ActiveView.FocusMap.Layer[0]).Workspace;

This works if the table of contents has no group layers.

Thanks


texts: Use Exiftools to add a tag with spaces and special characters to an image file

texts:
When I try to upload an image to Mapillary using it's web page based uploader and it's executable I get

EXIF Time issue.

When I compare a file that uploads correctly and the issue jpg file in EXIFTOOLS I see that it doesn't have Create Date, Date/Time Original and Modify Date. There are lots of other fields as well but these are the only ones dealing with time.

When I try to use

'.\exiftool(-k).exe' -"Date/Time Original"="2019:10:27 18:14:10.5" Photo_2019_Jul_12_13_38_40_019.jpg

It doesn't let me add the tag (I assume) because it has a '/' in it. It also doesn't allow "Create Date" for example.

How can I create a tag with these?

Some sample files are in https://drive.google.com/drive/folders/1QCQdSvdk0RygfCqqRWaOj-IzBbj17gax?usp=sharing

This is just to see if adding this will allow it to work in Mapillary but if anyone has an alternative solution to fix these hundreds of thousands of images that I have to upload it will be fine.


texts:
I have a QGIS project with 3 layers:

  OpenStreetMap layer
  Polygon representing a specified area on the map
  Lines representing streets in that area on the map

How do I convert this project to an OSM pbf file?
I'm trying to then pass that file to osm4routing  to process the street network and output a graph ready for routing.


texts: Calculate cells from which the feature would be visible above the skyline (horizon)

texts:
I'm trying to work out whether there is a way to do something similar to a ZTV/Visibility analysis but only include the cells from which the point (with 50m OFFSETA/height) would be visible above the horizon.

For example a windfarm sited at the top of a ridge may be said to have a large Zone of Theoretical Visibility in one direction, but if there is a mountain behind it then the visual impact of that turbine/farm would be lessened.

I would be grateful if anyone knows a way to do this. I can use any open source software or ArcGIS with Spatial Analyst and 3D licences.


texts: ArcGIS JS API Multiple Select In Where Clause of Query()

texts:
I have  a simple query like

query.where = "NAME = 'Renfrew'";

which works perfectly on QueryTask(). Now I would like to select array of points like

query.where = "NAME = 'Central Branch' AND NAME = 'Renfrew'";

but it is not working. Can you please let me know how to fix this?


texts: How do I transform an EPSG:3857 world map to an EPSG:4326 ECE world map?

texts:
This follows up a question posted over seven years ago. The results look very similar to the original inquiry without an answer or update I could test successfully.

I made a custom stylized EPSG:3857 map that matches up perfectly (with exception of a few differences in Antarctica) to an OpenStreetMap EPSG:3857 export:

My goal is to create an ECE image that lines up with NASA's Blue Marble (EPSG:4326 I think):

Image courtesy of NASA.

Following steps I found here, I can get close.

gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
-180 85.0511

Generates:

-20037508.3427892 20037471.2051371 0

Convert:

gdal_translate -of Gtiff -co "tfw=yes"  -a_ullr -20037508.3427892 20037471.2051371 20037508.3427892 -20037471.2051371 -a_srs "EPSG:3857"  "/espg-3857.tiff" "tfw.tiff"

gdalwarp -s_srs EPSG:3857 -t_srs EPSG:4326 -ts 1024 512 "tfw.tiff" "espg-4326.tiff"

Result:

This is not the desired result however. Note the closeness of Greenland to the top compared to the Blue Marble image. When overlaying the two, many differences can be seen:

What can I do to correctly transform from EPSG:3857 to EPSG:4326 for this world map?


texts: Updating/overwriting multiple layers within a feature layer in ArcGIS Online from a fgdb containing all the new feature classes

texts:
I am looking for a method to upload and replace multiple layers in ArcGIS Online that are stored within a (hosted) feature layer. I have the  feature classes to be used to update (and identical in name) stored in a fgdb on my local drive.
Does anyone have any ideas or tips to a script that are able to perform what I am looking for? I am fairly new to scripting, and have not been able to find other suggestions that updates multiple layers in a feature layer and deletes the old ones.


texts: Adding more than one legend in QGIS Print Composer?

texts:
I'm using QGIS 2.8.4 and I can't seem to figure out adding more than one legend to my Print Composer. Say I want a larger map and then an inset map -- showing perhaps a different classification/color scheme -- on the same page of my Print Composer. For whatever reason, despite checking the lock button(s), my "locked" legend will change according to whatever my original Q map shows. Therefore, if I want to add two legends, I'm really just adding the most-recent one and duplicating it.

I've tried every possibility of the locked buttons/cache/etc.


texts: Merging GDAL tiles and crop via bounding box?

texts:
I am using GDAL.

Assume I have four adjacent GDAL tiles in the formation of a square. I have a bounding box with a corner in each tile so that the rectangle formed spans all four tiles. My goal is to generate a DEM from those four adjacent GDAL files that is cropped to the bounding box. Also as an added bonus, assume these GDAL tiles are in folder of GDAL tiles. In this case, the DEM is a Geotile.

In essence I am asking how to perform three steps here:

Find all GDAL tiles covered by the bounding box.
Merge the tiles into one giant DEM
Crop the DEM to the bounding box and output as GEOTIFF.

The methods I know to do this would be really efficient (ex: open each GDAL sequentially and find all the ones that contain a corner point, merge tiles, then crop, surely GDAL has a more efficient means of doing this).


texts: How to know if a point is inside a rectangle on the map?

texts:
I have four points that represent a rectangle (a,b,c,d). I got these points from mobile GPS.

In my database I have this table of business locations:

bussId | longitude | latitude
-----------------------------
  ..      ..          ..

How can I know which businesses are inside the rectangle?

I want to implement this function:

function isInsideRectangle(Point buss,Point a,Point b,Point c,Point d)

where the Point contains longitude and latitude.


texts: How to transform JSON with field x,y to GDB in FME

texts:
I have JSON file with x,y field. I wonder how to transform geometry to GDB?

Example of JSON:

[
{
    "OBJECTID": "1",
    "x": "15",
    "y": "21"
},
{
    "OBJECTID": "2",
    "x": "16",
    "y": "21"
},
{
    "OBJECTID": "3",
    "x": "16",
    "y": "21"
}
 ]

I try use GeoJSON and JSON reader. JSON reader read attributes but I don't now how to insert x,y to shape field in GDB.

Now I create table with field contains WKT. How to convert table with field containing WKT to geometry layer?


texts: The first "n" points near a polygon? (spatialite or postgis query)

texts:
I've two spatial dataset, the fist one is a polygons layers (My_polygons) and the second one is a points layer (My_points).

I need to find, for each polygons, the "nearest" "N" points (N=10) using a spatial query in spatialite or postgis and store this information somewhere in my database for future queries.


texts: Adding multiple rasters produces unexpected results. I have set the data extent to the same in all rasters

texts:
I am trying to add three rasters together. I'm trying to figure out where I am going wrong so I am going to list out what I am doing:

I have Raster_1, Raster_2 and Raster_3.
Raster_1 and Raster_2 were set to the extent of Raster_3
I used the Map Algebra tool to write "Raster_1 + Raster_2 + Raster_3" = Raster_4
Raster_4 only has areas where all three previous rasters intersect.
I thought that by setting the extent of Raster_1 and Raster_2 to Raster_3 this wouldn't be an issue.
I went to check that I had done this correctly but the extents had gone back to their previous extents.
I reset Raster_1 and Raster_2 to the extent of Raster_3 but they are not the same. I have no idea why they would not be the same if that's how I am setting them. The setting I am using is:

  the rectangular extent of Raster_3

What am I doing wrong?


texts: How to change raster values of a specific area on the map

texts:
I have a raster file of a country surrounded by sea. I want the value, of the pixels that are above the sea (outside of the land), to be zero and everything else (pixels in the land) to have the initial value when i imported the raster file.
In the beginning every pixel has a random value (both, these on the land and on the sea). But I just want to zero all the values of the pixels on the sea


texts: How to remove cloudy pixel from MODIS NDVI (MOD13Q1)

texts: I am using MODIS NDVI product (MOD13Q1) for the vegetation changes study. I want to remove cloudy or bad pixel from MODIS NDVI data. I am not getting any suitable method or software or any program to do this task.


texts: Fastest way to perform focal operations using R

texts:
I know there are already some questions out there discussing this issue, for example:
How parallelize the extract function for raster files in R?
Increasing speed of crop, mask, & extract raster by many polygons in R?
But I'm not really sure if there is any "best practice" way to approach this. The setup of problem I describe is that I have one (or more, but the number of rasters does not matter too much I think) raster and a vector-object with many polygons. Now I want to get some aggregated value for each polygon.
What I found so for is that the normal raster::extract has the advantage of potentially returning a sp-object by setting sp=TRUE, however, it is raster slow.
The exact_extract-function seems to be way faster, but when using an aggregation function I can not append the aggregated value as a column at the vector (lets say sf-dataframe).
My question now would be if there is any way of making this even faster. Maybe using some kind of parallel approach or rasterizing the shape as proposed in the second link.


texts: Using GROUP BY linear spatial clusters with PostGIS

texts:
I have a series of points in a horizontal line, arranged into discrete spatial clusters by the underlying data pattern. Here's an example of what they look like:

How can I use PostGIS to group them spatially into the highlighted sets A, B and C?

This is a harder problem than it initially seemed to me. So far I've tried:

Using ST_X() to order them, then lag() to identify points between which there are gaps. Unfortunately this just groups the points into "gap" groups and "interior" groups.
Using convex hulls around arbitrary spatial thresholds. This both failed to separate all clusters and had the side effect of creating mid-cluster separations. (For future reference, this was a dumb idea.)
Using ST_GeoHash() to cluster them. This doesn't accurately catch the linear nature of the geometry.


texts: Map view does not render full SVG, but layers pane and print layout does?

texts:
In map view where I want to do all design and small changes the SVGs for "MELLANTID" does not render the added outline, only original size. The SVG for "SJUKVARD" only renders the left half of the SVG, only the drink symbol and not the medical.

I searched this site, Google and tried everything and still can't figure out why the SVGs render wrong in map view. They are correct in the layers pane and when exporting/print layout.


texts: Can anyone tell me  how to transform a shapefile into a WMS and then load it into mapbender?


texts: How can you use scikit-learn to do PCA on a 9 band raster image?

texts:
I would like to do a PCA on a 9band raster image to reduce the total number of inputs bands to 3. However, I'm not able to do it as I don't think I understand what sklearn is doing in its PCA.

Some code:

import skimage, rasterio, numpy, os
from rasterio.plot import reshape_as_raster, reshape_as_image

infile = r"C:/.../.../Sentinelimage9bands.tif"
im = rasterio.open(infile)
arr = im.read()
pca = sklearn.decomposition.pca(n_components = 3)
arrs = reshape_as_image(arr)
pca.fit(arrs)

Resulting error: "Found array with dim 3. Estimator expected <= 2"

How am I supposed to do the PCA if it can't accept more than 2 dimensions? Moreover, if I flatten the array into 1 dimension, how will I know the explained variance from each band?


texts: Writing multipart features in INTERLIS using FME

texts:
Trying to export a geopackage into an .xtf I run into the following problem:
unexpected geometry type COM.safe.fmeobjects.FMEMultiCurve
The .ili looks like this (simplified for posting here)
CLASS G (ABSTRACT) =
  Identifikator          : MANDATORY 1000 .. 9999;
  Name                   : MANDATORY TEXT*250;
UNIQUE Identifikator;
END G;

CLASS G_Linie EXTENDS G =
  Geometrie              : MANDATORY POLYLINE WITH (STRAIGHTS) VERTEX GeometryCHLV95_V1.Coord2;
END G;

GeometryCHLV95_V1 is a swiss base modul which you can find here:
It contains among other things:
Coord2 = COORD
  2460000.000 .. 2870000.000 [m] {CHLV95[1]},
  1045000.000 .. 1310000.000 [m] {CHLV95[2]},
  ROTATION 2 -> 1;

I assume the problems occur because the geopackage contains a multipart-Polyline and POLYLINE does not allow this.
I tried the following:
If I use an FME Deaggregator I run into problems because there is a UNIQUE identificator in the class.
I also tried the class Multiline
STRUCTURE MultiLine =
  Lines: BAG {1..*} OF LineStructure;
END MultiLine;

also part of the base model.
I applied it like this:
CLASS G (ABSTRACT) =
  Identifikator          : MANDATORY 1000 .. 9999;
  Name                   : MANDATORY TEXT*250;
UNIQUE Identifikator;
END G;

CLASS G_Linie EXTENDS G =
  Geometrie              : GeometryCHLV95_V1.MultiLine;
END Gefahrenquelle_Linie;

unlike the first 2 tries, this translation runs without errors, but the resulting .xtf file does not contain any Geometrie attribute on the written instances of G.
How do I get mulitpart features in my class?


texts: Finding centroid of point pairs based on common ID using ArcGIS Desktop?

texts:
I want to take a series of point pairs within a single shapefile and find a single centroid point for each pair.

For example in the image below, I would like there to be one point "31", one point "32", etc. without having to create a new point one by one.


texts: I'm wanting to create a service area within network analyst and want to input polygons as the facilities. However, I'm having trouble with this! I can input point features into the facilities, but am unable to do so with my polygon features. Can you only use point geometry as facilities within a service area analysis?


texts: Trying to use Symbols to represent numerous fields

texts:
I am using QGIS 2.12 in an Archaeological setting. I have populated a layer with features that represent burials and I have the intent for symbols to appear in the burial polygons to represent the type of artifacts(shell beads, types of pendants, etc.) found in them. I have 6 different items I want to represent on these. Some burials have 1 of these items, some burials have all 6.

How do I go about setting it up so that the fields can all be represented? Can I use rule based to do it?


texts: Looking for political maps (administrative division)

texts:
I'm looking for an interactive map that shows ONLY administrative divisions of a country. I need to personalize it so when the user clicks a region, the map shows some information about it. If it comes in colors and animations, it would be better. Is there a way I can achieve this? I've seen around and I found very good maps but maybe they show additional geographic information that is really good but might not be relevant for my web application (PHP).
One more thing: I'm avoiding the use of Flash and Java Applets.
Thanks.


texts:
I'm a new member and a beginner with Leaflet and JavaScript.

I try to display shapefile data on my webmap with Leaflet. But it doesn't work.
All the data (.shp, index.html... ) in WGS84 are on my computer, in the same directory.

My JavaScript code is:

var shapeLayer = new L.Shapefile("data/QualiAgro.zip");

I get the following Console message:

  Cross-Origin Request Blocked: The Same Origin Policy disallows reading
  the remote resource at https://some-url-here. (Reason: additional
  information here).

If I understood this is a CORS problem (url). But I don't know how to do?
Does any one have any idea?

I also try the shapefile data with: http://leaflet.calvinmetcalf.com/#3/32.69/10.55
There's no display.

Thank you Vince for your answer.

I have installed "Wampserver2" on my computer to run my html from a local server.
I changed the JavaScript code :

    var shapeLayer = new L.Shapefile("http//localhost:80/data/QualiAgro.zip");

But I still have the same problem with no message Console for the mistake.


texts: How to import multiple shapefiles each with different source projections to postgis

texts:
I have a couple thousand .zip files each with a .shp, .dbf, .shx, and .prj inside them that I want to load into postgis. The problem is, they have many different source projections (each defined in the .prj file).

Is there any way to load them all at once without having to manually set the projection for each file?


texts: OpenLayers setMap function is removing controls I have added

texts:
I'm trying to extent Ol.Control with my own controls, like this:

class ZoomHistory extends ol.control.Control {
    constructor(opt_options) {
        var options = opt_options || {};
        var controlDiv = document.createElement('div');
        controlDiv.className = options.className || 'btn-group mr-2';

        var backButton = document.createElement('button');
        backButton.className = options.backClassName || 'm-btn btn btn-secondary';
        backButton.textContent = options.backLabel || '';
        backButton.title = options.backTipLabel || 'Previous view';
        backButton.addEventListener('click', function (evt) {
            var history = _this.get('history');
            var index = _this.get('index');
            if (index > 0) {
                _this.setProperties({
                    shouldSave: false,
                    index: index - 1
                });
                _this.getMap().getView().setProperties(history[index - 1]);
            }
        });
        backButton.disabled = true;

        var img = document.createElement('img');
        img.src = '../images/zoom-last.png';
        img.alt = "Previous view";
        backButton.appendChild(img);

        controlDiv.appendChild(backButton);

        super({
            element: controlDiv,
            target: options.target
        });
        var _this = this;

        _this.setProperties({
            history: [],
            index: -1,
            maxSize: options.maxSize || 50,
            eventId: null,
            shouldSave: true
        });

        _this.on('change:index', function () {
            if (this.get('index') === 0) {
                backButton.disabled = true;
            } else {
                backButton.disabled = false;
            }
            if (this.get('history').length - 1 === this.get('index')) {
                nextButton.disabled = true;
            } else {
                nextButton.disabled = false;
            }
        });

    }
    setMap(map) {
        if (map === null) {
            ol.Observable.unByKey(this.get('eventId'));
        } else {
            this.set('eventId', map.on('moveend', function (evt) {
                if (this.get('shouldSave')) {
                    var history = this.get('history');
                    var index = this.get('index');
                    history.splice(index + 1, history.length - index - 1);
                    if (history.length === this.get('maxSize')) {
                        history.splice(0, 1);
                    } else {
                        index += 1;
                    }
                    history.push(map.getView().getProperties());
                    this.set('index', index);
                } else {
                    this.set('shouldSave', true);
                }
            }, this));
        }
    }
}

When I add this control in the map, it is not adding the control(backButton)

var map = new ol.Map({
    target: 'm_gmap_1',
    controls: [
        new GSPrint({ target: 'map_controls' }),
        new ZoomHistory({ target: 'map_controls' })
    ],

but when I comment setMap(map) function in  ZoomHistory class, the controls(backbutton) is getting added inside map_controls div, but button click event is not firing.


texts: How to add multiple layers to OpenLayers from a single Geoserver WMS?

texts:
I'm trying to add multiple layers from a single Geoserver source to my Openlayers map. I can get each layer to appear by changing a single line of code. but when I attempt to have both layers on at once, I get a blank map. How do I add both? One is a raster, the other is a vector.

I've attempted via having the additional layer in a new variable and using map.addlayer. This didn't break the map, but it still didn't load the additional layer. I tried by just adding the additional layer into the layer variable, broke the map. I also tried by adding it to the params. So it was

params:{'LAYERS': 'Test:TestImagery', 'Test:Vector'},

This broke the map. How do I make this work?

Code:

<script>
  var layers = [
    new ol.layer.Tile({
      source: new ol.source.OSM()
    }),

    new ol.layer.Image({
      extent: [-20000000, 1000000, -5000000, 9000000],
      source: new ol.source.ImageWMS({
        url: 'http://localhost/geoserver/Test/wms',
        params: {'LAYERS': 'Test:TestImagery'},
        ratio: 1,
        serverType: 'geoserver'
      })
    })

  ];
  var map = new ol.Map({
    layers: layers,
    target: 'map',
    view: new ol.View({
      center: [0,0],
      zoom: 14
    })
  });

</script>


texts:
I am using gdalwarp to reproject a .tif file to the same crs as a shapefile I want to clip the raster by later on.
When doing gdalwarp using

NDVI = gdal.Open("NDVI_2019.tif")
NDVI_rep = gdal.Warp("NDVI_repr.tif",NDVI,dstSRS = crs)

where crs is the crs of my Shapefile, the following output (right) is created out of the left image.

As you can see there is this "fill" to get a rectangle, where every value is just 0. Is there a way to change the values around my actual tif to something like -1000 to filter it?

Thanks in advance,
Felix


texts: How to convert a spatial Oracle LOB object type to a Geodataframe listlike geometry object column

texts:
I have in my dataframe one column with cx_Oracle.LOB object [2D polygon type], made with a SDO_UTIL.TO_WKTGEOMETRY SQL function and I want to create a geodataframe, and set my geometry column the geometry returned from the SDO_UTIL.TO_WKTGEOMETRY function to a valid geodataframe listlike geometry.

The Geodataframe function return error when i set the geometry column with the WKTGEOMETRY made with SDO_UTIL.TO_WKTGEOMETRY function. I work with PyCharm application.

my dataframe:

dfparc = pd.DataFrame(dfparc.values, columns=['PAR_ID', 'PAR_NUM', 'geometry', 'AREA_HA', 'DCF'])

print(type(dfparc.geometry[0]))

class 'cx_Oracle.LOB'>

print(dfparc.geometry[0])

POLYGON ((225215.22 411994.91, 225214.02 411912.51, 225339.62 411916.91, 225339.62 411995.31, 225279.22 411995.31, 225215.22 411994.91))

My Geodataframe:

geoparc = gp.GeoDataFrame(dfparc, crs=102160, geometry=dfparc.geometry)

returned errors:

...Python35-32\lib\site-packages\geopandas\geodataframe.py", line 50, in __init__
     self.set_geometry(geometry, inplace=True)
...Python35-32\lib\site-packages\geopandas\geodataframe.py", line 149, in set_geometry
    raise TypeError("Input geometry column must contain valid geometry objects.") TypeError: Input geometry column must contain valid geometry objects.

I suppose the problem could be in the type of object returned from SDO_UTIL.TO_WKTGEOMETRY function that is not available for geodataframe geometry data type.


texts: Clipping line layer based on polygons with QGIS?

texts:
I have a line layer and a polygon layer. Both files have the same SRS. I want in a new shape file only the lines covered by the polygons.

I tried in QGIS the Vector -> Geoprocessing -> Clip tool but I only get an empty shapefile.

How do I do this with QGIS?


texts: Calculate field values using arcpy.da.UpdateCursor

texts:
I am trying to create and update a field with a count of line features (tLayer) within a distance of point features (sLayer). I am attempting to use a combination of AddField_management, arcpy.da.SearchCursor, SelectLayerByLocation_management, arcpy.GetCount_management, and arcpy.da.UpdateCursor.

The code I have for this is currently updating all records for the Line_Count field with the count of the point features (i.e. 2) for only the last record. Though, a print statement following the GetCount line will return the line count for all of the point features (with a few unessential iterations). I believe it is only printing the last row in 'result' for the Line_Count update statement. Is there an efficient way for result (i.e. result=int(arcpy.GetCount_management("tLayer_lyr").getOutput(0))) to be saved as a list or an array that can be written to the Line_Count field? ...Or would it be best to use a dictionary?

What do I need to do to appropriately update the Line_Count field for all of the records? Also, this process will be applied to a large dataset and will be extended to include 'where clauses'. Are there any suggestions as to how to make this as efficient as possible?

import arcpy
from arcpy import env
arcpy.env.OverwriteOutput = True

defaultGdbPath = 'C:\Topo_Check_Tess_V5.gdb'

sLayer='C:\Topo_Check_Tess_V5.gdb\Stations'
tLayer='C:\Topo_Check_Tess_V5.gdb\Lines'
#ppLayer='C:\Topo_Check_Tess_V5.gdb\Plants'

arcpy.AddField_management(sLayer, "Line_Count", "SHORT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

TLineCountField = "Line_Count"

arcpy.MakeFeatureLayer_management(tLayer, "tLayer_lyr")
with arcpy.da.UpdateCursor (sLayer, TLineCountField) as LineCountList:
       for s_row in LineCountList:
            with arcpy.da.SearchCursor(sLayer,["OBJECTID", "SHAPE@"]) as s_list:
                    for station in s_list:
                        arcpy.SelectLayerByLocation_management("tLayer_lyr", 'WITHIN_A_DISTANCE', station[1],.000002, "NEW_SELECTION")
                        result=int(arcpy.GetCount_management("tLayer_lyr").getOutput(0))
                        print result
                        LineCountList.updateRow(s_row)
                    del station
            del s_row

print "Done"


texts: Exporting EVI image on Google Earth Engine with visual parameters between -1 and 1?

texts:
I am having an issue with exporting an EVI image with the visual parameters between -1 and 1. On Google Earth Engine (GEE), using the inspector tool, the values are between -1 and 1, but when I export and upload onto Arcmap, the values exceed that threshold to something like -300 to 800. This isn't an issue with the NDVI.
var image = ee.Image('LANDSAT/LC08/C01/T2_TOA/LC08_097055_20180301');

var evi = image.expression(
    '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
      'NIR': image.select('B5'),
      'RED': image.select('B4'),
      'BLUE': image.select('B2')
});

Map.centerObject(image, 9);

var eviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};

Map.addLayer(evi, eviParams, 'EVI image');

Export.image.toDrive({

  image: evi.visualize(eviParams), //I have tried solely with **evi** and that doesn't work either

  description: 'EVIex',

  scale: 30,

  region: geometry

});


texts: Stop arcpy.CopyFeatures_management from adding copied feature classes to TOC

texts:
I've written a script that copies several feature classes to a new file gdb. Unfortunately CopyFeatures_management automatically adds the new feature classes to the TOC.
Is there a way to suppress this behaviour?
import arcpy, os

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
layers = arcpy.mapping.ListLayers(mxd,"",df)

from arcpy import env
env.overwriteOutput = True
outgdb = r"C:\\temp\\test.gdb"
for lyr in layers:
    if lyr.isFeatureLayer:
        outfc = os.path.join(outgdb,os.path.split(lyr.dataSource)[-1])
        arcpy.CopyFeatures_management(lyr.longName, outfc)


texts: How to programmatically associate a layer with a python init method from external file?

texts:
I know it is possible to associate a layer with a custom UI for editing feature attributes using the QGIS interface by using python:

vectorLayer1.setEditForm('/home/me/uifile.ui')

And I found this to associate the vector layer with python init method:

vectorLayer1.setEditFormInit('moduleName.functionName')

But this only allowed when the module (python file) located in the same directory as my project. Is there any ways to load from external file?


texts:
i get the following Error when running the code: ImageCollection (Error) User memory limit exceeded.

I want to calculate the difference between two images in an image  Collection.
Can anybody help me?
Here is my code:

// Load Sentinel-2 TOA reflectance data.
var Sentinel2 = ee.ImageCollection('COPERNICUS/S2')
              .filterDate('2018-01-01', '2018-12-31')
              // Pre-filter to get less cloudy granules.
              .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
              .map(maskS2clouds)
              .filterBounds(Greece);

function addNBR(image) {
// Compute the Normalized Difference Vegetation Index (NDVI).
var nir = image.select('B8A');
var swir = image.select('B12');
var nbr = nir.subtract(swir).divide(nir.add(swir)).rename('NBR');
return image.addBands(nbr)
}
var Sentinel2_collection = Sentinel2.map(addNBR)

print ('Sentinel2_collection: ', Sentinel2_collection)

var NBR = Sentinel2_collection.select('NBR').sort('system:time_start');

print ('NBR is ', NBR)

var NBRlist = ee.List(NBR.aggregate_array('system:time_start'))
print ('NBRlist ', NBRlist)

var NBRdiff = ee.ImageCollection.fromImages(NBRlist.map(function(date){

  // get the image corresponding to the date
  var currentImage = NBR.filter(ee.Filter.eq('system:time_start',    date)).first();

 // Now we have to get the 'previous' image in the collection
  var indexCurrent = NBRlist.indexOf(date);
  var indexPrevious = indexCurrent.subtract(1);
  var datePrevious = NBRlist.get(indexPrevious);
  var previousImage = NBR.filter(ee.Filter.eq('system:time_start',    datePrevious)).first();

  // Subtract the current image from the previous
  var diffImage = currentImage.subtract(previousImage).select(['NBR'], ['NBRdiff']);

 // return the image
  return currentImage.addBands(diffImage)
        .set('system:index_previous', previousImage.get('system:index'));
}));

// print and add to the map
print('Difference Image Collection', NBRdiff);


texts: FME Desktop Distorts shapefile to TDRBM(.bas) translation

texts:
I am trying to convert esri shapefile into Tobin TDRBM II Data Distribution Format (.bas) but this distorts my data and create files with extension .bas.csv instead of .bas.I removed .csv and tried to import in global mapper.My shapefile is in WGS84. I tried both specifying and not specifying the coordinate system in the FME but all went in vain.

Before translation-

After translation-

Even i lost my attributes in the shapefile after translation!

My FME details:

Edition: FME Desktop Professional Edition (floating)
Version: FME(R) 2016.1.0.0 (20160419 - Build 16492 - WIN64)
Locale: en_US
Codepage: 1252 (ANSI - Latin I)


texts: Georeference map with just one known point

texts:
I have an old scanned topographic map I have to georeference. Unfortunately there is no coordinate system stated and it has just on reference grid point on it. Besides that it states in the down left corner x and y in 1000 steps ranging from 1000 to 3000. It also states the scale 1:10000.

What would be the best way to georeference this kind of map? Would the scale be of any use? It is a mountainous area too, so the surface changed a lot in the meanwhile.

Edit: Here is a small part of the map for better illustration. It could be, that they set up their own local reference system, so would it be best to georeference using features on the map?


texts:
I am attempting to use GDAL/OGR to read from an SDE enabled database (MS SQL). I have tried the steps located here but do not have access to the database server. I am also using a GDAL installation for python downloaded and installed from GISInternals. (gdal-202-1500-core.msi and GDAL-2.2.1.win32-py2.7.msi) so my setup is a little different.

What I have tried

I downloaded the archived ms4w_3.0.6 and copied the ogr_SDE.dll and gdal_SDE.dll to C:\Program Files (x86)\GDAL\gdalplugins
Downloaded and installed the ArcSDE binaries from the MyEsri site. Verified  the SDEHOME path was set correctly.
Add C:\Program Files\ArcGIS\ArcSDE\sqlexe\bin;C:\ms4w_3.0.6\ms4w\tools\gdal-ogrto path.

However, when I check for available formats, I see this:

C:\Program Files (x86)\GDAL\gdalplugins>ogrinfo --formats
ERROR 1: Can't load requested DLL: C:\Program Files (x86)\GDAL\gdalplugins\ogr_S
DE.dll
126: The specified module could not be found.

ERROR 1: Can't load requested DLL: C:\Program Files (x86)\GDAL\gdalplugins\ogr_S
DE.dll
126: The specified module could not be found.

Supported Formats:
  -> "ESRI Shapefile" (read/write)
  -> "MapInfo File" (read/write)
  -> "UK .NTF" (readonly)
  -> "SDTS" (readonly)
  -> "TIGER" (read/write)
  -> "S57" (read/write)
  -> "DGN" (read/write)
  -> "VRT" (readonly)
  -> "REC" (readonly)
  -> "Memory" (read/write)
  -> "BNA" (read/write)
  -> "CSV" (read/write)
  -> "NAS" (readonly)
  -> "GML" (read/write)
  -> "GPX" (read/write)
  -> "LIBKML" (read/write)
  -> "KML" (read/write)
  -> "GeoJSON" (read/write)
  -> "Interlis 1" (read/write)
  -> "Interlis 2" (read/write)
  -> "GMT" (read/write)
  -> "GPKG" (read/write)
  -> "SQLite" (read/write)
  -> "ODBC" (read/write)
  -> "WAsP" (read/write)
  -> "PGeo" (readonly)
  -> "MSSQLSpatial" (read/write)
  -> "PostgreSQL" (read/write)
  -> "MySQL" (read/write)
  -> "PCIDSK" (read/write)
  -> "OpenFileGDB" (readonly)
  -> "XPlane" (readonly)
  -> "AVCBin" (readonly)
  -> "AVCE00" (readonly)
  -> "DXF" (read/write)
  -> "Geoconcept" (read/write)
  -> "GeoRSS" (read/write)
  -> "GPSTrackMaker" (read/write)
  -> "VFK" (readonly)
  -> "PGDump" (read/write)
  -> "OSM" (readonly)
  -> "GPSBabel" (read/write)
  -> "SUA" (readonly)
  -> "OpenAir" (readonly)
  -> "PDS" (readonly)
  -> "WFS" (readonly)
  -> "HTF" (readonly)
  -> "AeronavFAA" (readonly)
  -> "Geomedia" (readonly)
  -> "EDIGEO" (readonly)
  -> "GFT" (read/write)
  -> "GME" (read/write)
  -> "SVG" (readonly)
  -> "CouchDB" (read/write)
  -> "Idrisi" (readonly)
  -> "ARCGEN" (readonly)
  -> "SEGUKOOA" (readonly)
  -> "SEGY" (readonly)
  -> "ODS" (read/write)
  -> "XLSX" (read/write)
  -> "ElasticSearch" (read/write)
  -> "PDF" (read/write)
  -> "Walk" (readonly)
  -> "CartoDB" (readonly)
  -> "SXF" (readonly)

C:\Program Files (x86)\GDAL\gdalplugins>

I am simply looking to read data from the SDE database with GDAL/OGR.

Looks like I am trying to connect to ArcSDE version 10.2.1. Not sure if its even possible to connect to a 10+ SDE with GDAL.


texts: How to label pie charts in ArcGIS for Desktop?

texts: How can I add to pie charts in ArcGIS (10.1) labels of the values the charts are created with?


texts:
I was wondering if there was any python (or R) packages out there that can do spatial clustering on ranked data?

For each US state I spatially varying value for education rankings. For example, for Virginia the ranking may be 50 compared to a ranked value of 10 for the state of Missouri. These are not categorical (e.g. A vs B) or ratio (e.g. 1.5), but ranked in the sense that a value of 50 for Virginia does not mean that this state is 5 times higher in educational value that Missouri. Its synonymous to ranked outcomes at the end of a marathon (e.g. 1st, 2nd,....places). Most of the stuff I've been finding online relate to the relationship between records across different attribute themes in tabular data (e.g. Gower distance). I'm looking for a spatially clustering method such as Moran's I where the attribute vector being analyzed is an education rank value for each state polygon.


texts:
I'm frequently working with the imagery basemaps provided by Esri.

I find that the labels (cities, districts etc.) are way too small.

Is there any option to increase the font size?


texts: How do I use the True Type Font that has a space in the name in an SLD?

texts:
Using an answer from another question, I learned how to use a Wingdings character in my point SLD.  Now, I'd like to use a character from "Wingdings 3".

Here's the XML from the answer above that works great with Wingdings.  (I switched the character to 116 as I want a triangle.

<Graphic>
  <Mark>
    <WellKnownName>ttf://Wingdings#116</WellKnownName>
    <Fill>
      <CssParameter name="fill">#000000</CssParameter>
    </Fill>
  </Mark>
  <Opacity>1</Opacity>
  <Size>18</Size>
</Graphic>

Some of the combinations I've tried.

ttf://Wingdings3#116
ttf://Wingdings 3#116
ttf://Wingdings$203#116
ttf://Wingdings+3#116
ttf://Wingdings_3#116
ttf://WINGDNG3#116

What combination am I missing?


texts: Adding two rasters results in negative infinity and positive infinity

texts:
I'm trying to create a hydro-flattened DEM according to USGS specifications. I have two rasters, a DEM of my study area and a raster that only contains the water bodies that I have flattened so the pixels represent a flat surface for each of the lakes.

I've created a polygon mask that I've used to remove the lakes from the original DEM and used the saga:invertnodata tool to create an opposite mask that I used to generate the second raster containing only the flattened lakes.

When I try to combine the two rasters, my resulting image is blank with the only values recorded in the legend as 1.79769e+308 and -1.79769e+308 (for some reason -infinity is the upper bound but that may be trivial at this point).

I cannot for the life of me figure out why adding the two rasters (using simple addition in raster calculator) results in -infinity and +infinity. I suspect is has to do the masked values in each raster have a value of 'no data' but I've tried converting each raster's 'no data' values through gdal warp, translate, and a couple other tools that are supposed to convert no data values to a specified value.

Has anyone ever encountered this, or have another method to handle and convert 'no data' to 0?


texts: Renaming the result of a QGIS Processing Algorithm

texts:
When calling a QgsAlgorithm from within another, the result will be named after this algorithm output regardless of what is setup as the output of the main one.

For instance, taking the code from the answer to this question, although the output is set with the name 'OUTPUT' and the returned dictionary is the following:

{'OUTPUT': 'clipped_xxxxxxxxxxxxxxxxxxxxx'}

the output layer is actually named 'clipped' and I cannot seem to be able to change this.

I've tried the setName() method on the QgsVectorLayer retrieved with QgsProcessingUtils.mapLayerFromString() but to no avail.


texts: Heatmap corrected for population density

texts:
I would like to make something like a heatmap, but have it corrected for another factor, say population density.

For example, lets say I want to create a heatmap that displays whether or not somebody has a dog.  People in low population density areas are more likely to have a dog, but a basic heat map will make higher population areas look like they have more dogs just because there is a lot more people.

I would like to somehow make the heatmap take this factor into account so it shows a higher value for the # dogs by population density


texts:
I'm very new to GIS and for learning I started with building a simple web app with GeoDjango. I am using dango-leaflet. Since I have a very limited knowledge I am facing many challenges. I'm trying to put a specifc marker on my map as a test
models.py
from django.db import models
from django.contrib.gis.geos import point
from django.contrib.gis.db import models as gismodels

class Points(gismodels.Model):
      data_lat = 44.915223
     data_long = -93.209741

     @property

     def geom(self):
         return point(self.data_lat, self.data_long, srid=4326)
     def __unicode__(self):
        return self.title

about.html
<body>
      <h1> Map  </h1>
     {% leaflet_map "main" callback="main_map_init" %}

      <script type="text/javascript">

          function main_map_init (map, options){

                var dataurl = '{% url "data" %}';

               $.getJSON(dataurl, function(data) {

                      L.geoJson(data).addTo(map);
               })
          }
      </script>

</body>

I'm having issue with data not being passed on to map. It would be nice if someone could point me to the right direction.


texts: Can I save an ArcGIS geospatial PDF with layers turned off, which can then be turned on individually in a pdf viewer?

texts:
I am working on an mxd that has many individual point layers (displaying species locations). To make life easy for my stakeholders, I'd like to export to a geospatial pdf so that all of these point layers are turned off when they open the pdf. They can then select which species they wish to view rather than having to turn off ALL species layers to start with.

Currently I am only able to save to geospatial pdf with all layers turned on. When I turn off all layers in my toc in arc, then save to geospatial pdf, those layers aren't saved to the pdf.

Any ideas?


texts: Can percent overlap of polygon layers be calculated in geographic coordinate system?

texts:
I have two global polygon shapefiles in WGS84, and I would like to generate spatial stats on the % area of one shapefile that is overlapping the other. I understand that geometric calculations such as area require a projected coordinate system, but I am wondering if that is necessary in order to calculate % overlap.

I would like to avoid projecting, as these are large datasets. However, in the event that I do have to project, are there recommended projections for running geometric calculations on global data? I've only ever worked with local / regional projections. I am using ArcMap 10.2.


texts: Web Mapping Client with front-end query function

texts:
I'm searching for a cloud based Web Map Client that offers a front-end query tool so that my client can query against the data and visualize the results on the web map.  I'm trying to keep costs down and looking at Open Source and lower tiered subscriptions. I have no experience with development and coding. Research so far:

Mango Map - query tool available for Professional tier
CartoDB - SQL query on back end, not sure if it's available on the front-end
MapGuide - has a query tool, but uncertain how to setup a server to host it

ArcGIS Explorer Desktop is a wonderful solution and an excellent example of what I'm trying to achieve, but the client is running a Mac OS and therefore cannot use the application without a Business account.


texts:
I'm trying to join one column into another.  one is "acct id" and the other is "taxid number".  Acct is from one county whereas the tax is from another, so there are blanks in the number of parcels/rows.  I did it before last week and I had it written down how to do this.

What I need to do is move 'acct id' into the 'taxid' column. I remember having " " as part of the vb script.  What are the steps involved as I keep getting a blank field when I attempt it.


texts: Combining ListDatasets and ListFields to return geodatabase contents to CSV--one script prints attributes and another returns feature class counts

texts:
New to Python. I'm trying to combine two scripts I found online that use the ListDatasets and ListFields tools. One loops through all the features of a geodatabase and writes the name of the feature dataset, feature class, and the number of features to a .csv file. To that one I would like to add the other that prints to the screen information about the attributes--name, type, and length. I would like all of this printed to a single .csv file like this:

The ListDatasets Script:

import arcpy

arcpy.env.workspace = "C:\Users\Redfish\Test.sde\Test.DBO.Structures\Test.DBO.HospitalBldg"

file = open("C:\Users\Redfish\Desktop\DatabaseDictionary.csv", "w")
file.write('Dataset, LayerName, NoOfRecords' + '\n')

for fds in arcpy.ListDatasets('*','feature') + ['']:
                for fc in arcpy.ListFeatureClasses('*','',fds):
                                result = int(arcpy.GetCount_management(fc).getOutput(0))
                                file.write(fds + ', ' + fc + ', ' + str(result) + '\n')

file.close()

The ListFields Script

import arcpy

arcpy.env.workspace = "C:\Users\Redfish\Test.sde\Test.DBO.Structures\Test.DBO.HospitalBldg"

fields = arcpy.ListFields("C:\Users\Redfish\Test.sde\Test.DBO.Structures\Test.DBO.HospitalBldg")
for field in fields:
    print("{0}\{1}\{2}"
          .format(field.name, field.type, field.length))


texts:
I created a WMS in a QGIS project and serve it with QGIS Server.
When I click on a feature on the WMS map, no attribute information is displayed.
I guess, the attributes are not exposed to the WMS but I do not know how to do this.
What settings are needed to expose attributes in QGIS to a WMS?


texts:
I'm trying to connect to a wfs using GDAL, from reading some of the (fairly sparse) documentation on this i have found that this is possible if GDAL is compiled with the correct drivers (Curl, Xerces).

I have compiled GDAL against both of these and have a working WFS Service to connect to.

From reading the documentation it says this is how the connection should be opened.

GDALDataset  *poDataset;

GDALAllRegister();

poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );

I have done this and replaced pszFilename with a URL to my WFS which didnt work and gave the error (and i have no idea what file.dat is referring to as all im doing here is entering a URL that should return some XML):

ERROR 4: `/vsimem/http_1/file.dat' not recognised as a supported file format.
ERROR 4: Failed to create temporary file:/tmp/file.dat

alternatively it says that you could use an XML File that contains the URL like this:

<OGRWFSDataSource>
    <URL>http://path/to/WFS/service[?OPTIONAL_PARAMETER1=VALUE[&amp;OPTIONNAL_PARAMETER2=VALUE]]</URL>
</OGRWFSDataSource>

I have made this XML file and put in my URL for the WFS im trying to connect to and i get this error even though the documentation says that reading the url from XML is supported:

ERROR 4: `:/Resources/WFS.xml' does not exist in the file system,
and is not recognised as a supported dataset name.

Im wondering if anyone knows how to connect to a WFS using GDAL because as far as documentation goes this is about all its says in the way of how to connect and goes into no detail about supported file formats or how to connect to any of these.

Any help would be appreciated (oh and btw my WFS is working because if i type the URL into my web browser i get the XML (or is it GML?) in the web browser for the request i send (a GetCapabilities one))


texts: Shapefile attribute table to Datagridview with C#

texts: Can you please let me know how to display a shapefile attribute table in C# DataGridView?


texts: GeoServer intersect filter error with SQL view publish

texts:
I have created a GeoServer service with a SQL view. The service is serving rivers. I want to select rivers in a polygon. So I am trying GeoServer demo request tool.

<wfs:GetFeature service="WFS" version="1.0.0"
  outputFormat="GML2"
  xmlns:topp="http://www.openplans.org/topp"
  xmlns:wfs="http://www.opengis.net/wfs"
  xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/wfs
    http://schemas.opengis.net/wfs/1.0.0/WFS-basic.xsd">
  <wfs:Query typeName="topp:rivers">
    <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
      <ogc:Intersects>
        <ogc:PropertyName>geom</ogc:PropertyName>
        <gml:Polygon xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4321">
          <gml:exterior>
            <gml:LinearRing>
              <gml:posList>78.57 58.32 78.48 24.7 96.89 81.23 51.61 62.41 08.42</gml:posList>
            </gml:LinearRing>
          </gml:exterior>
        </gml:Polygon>
      </ogc:Intersects>
   </ogc:Filter>
</wfs:GetFeature>

This request gives an error:

java.lang.RuntimeException: java.io.IOException java.io.IOException null ERROR: Operation on mixed SRID geometries

If my request goes to a service that not created with aSQL view, it works. How can I solve this problem?


texts:
I want to save the geometry data after I modify it, so I need an event that fire after I finish to modify the feature.
I use the following from this post:
listener modify interaction

    var selected_features = select.getFeatures();
// when a feature is selected...
selected_features.on('add', function (event) {
    // get the feature
    var feature = event.element;
    feature.on("change", function (event) {
            logStatus(feature.getId() + " has been changed!\n");
    });
});

But this listener fires too many times in a single modification, and if I am using 'once' instead of 'on' its fires only the first little change and no more, so its not saving the real change.
I need something like "changeend" or "modifyend" that will fire when I will finish to modify the feature, what can i use?


texts: How do I use code (preferably c#) to clip all points within a shapefile that are not contained within a polygon?

texts:
So the task is simple: I want to take a point shapefile as input and remove all of the points that are not contained within a polygon.

I'm using DotSpatial to do this, but the catch is that there are millions of points, and it's just... hilariously slow.

I'm using some code that looks like this:

for (int i = 0; i < points.ShapeIndices.Count - 1; i--)
{
    if (polygons.ShapeIndices[activeRow].Intersects(points.ShapeIndices[i]))
    {
        allPoints.Add(points.Features[i]);
    }
}

This takes minutes to do 80,000 points. There's obviously something janky with their implementation but I don't have time to run it down.

So my question is, if you had to get this done and introduce minimal dependencies, how would you do it? I could install PostGIS and load the shapefiles into PostGIS to do the contains operation, but I am loath to introduce a database into this project.


texts:
I have two sets of data, one is in OS coords, no problem there, but the other was set out in about 50 off 10x10 metre grids aligned to a 'feature' on the ground. I do have the OS coordinates of their datum, but that's it.

I'd appreciate some help on how to try and convert these grid data points to OS grid coords.

I am using QGIS 1.8, but have access to a CAD system if that is easier.

The four reference OS coords are:
* 552369 265951
* 552363 266159
* 552470 266068
* 552625 265929

and the grid for the data are in the form, where x and y are in metres from the datum line
* X Y   No.
* 600   20  2
* 600   30  6
* 600   40  0
* 610   10  4
* 610   20  4


texts: Layer ordering of ST_ASMVT tiles in Leaflet

texts:
I am generating MVT tiles for multiple layers in Django using the following query:
WITH
    bounds AS (
            SELECT ST_Segmentize(ST_MakeEnvelope(8247861.100083642, 3551570.082242422, 8252753.069893893, 3556462.052052673, 3857),1222.9924525627866) AS geom,
                            ST_Segmentize(ST_MakeEnvelope(8247861.100083642, 3551570.082242422, 8252753.069893893, 3556462.052052673, 3857),1222.9924525627866)::box2d AS b2d
    ),
    f_282 AS (
            SELECT ST_Simplify(ST_AsMVTGeom(ST_Transform(t.geom, 3857), bounds.b2d, 4096, 128, true), 25) AS geom,
            t.id
            FROM g_feat_poly t, bounds
            WHERE ST_Intersects(t.geom, ST_Transform(bounds.geom, 4326))
            AND t.shp_id = 282 AND t.deleted = FALSE
            UNION
            SELECT ST_Simplify(ST_AsMVTGeom(ST_Transform(t.m_geom, 3857), bounds.b2d, 4096, 128, true), 25) AS m_geom,
            t.id
            FROM g_feat_poly t, bounds
            WHERE ST_Intersects(t.m_geom, ST_Transform(bounds.geom, 4326))
            AND t.shp_id = 282 AND t.deleted = FALSE
    ),
    f_259 AS (
            .........
    ),
    f_258 AS (
            .........
    ),
    f_255 AS (
            .........
    ),
    f_257 AS (
            .........
    ),
    sa AS (
            .........
    )

    SELECT ST_AsMVT(f_282.*, '282') FROM f_282
        UNION
    SELECT ST_AsMVT(f_259.*, '259') FROM f_259
        UNION
    SELECT ST_AsMVT(f_258.*, '258') FROM f_258
        UNION
    SELECT ST_AsMVT(f_255.*, '255') FROM f_255
        UNION
    SELECT ST_AsMVT(f_257.*, '257') FROM f_257
        UNION
    SELECT ST_AsMVT(sa.*, 'sa') FROM sa

and displaying it in Leaflet using vectorgrids with following code:
var url = '..../api/v1/tiles/6/13/5782/3369/?shp=282-POLYGON,259-POLYGON,258-POLYGON,255-POLYGON,257-POLYGON,7-sa'
var featTilesLayer = L.vectorGrid.protobuf(url, {
                interactive: true,
                vectorTileLayerStyles: shpStyle,
                getFeatureId: function(f) {
                    return f.properties.id;
                }

But i am unable to change the order of layers as all the layers are in single vector layer. I tried changing the order in URL and in SQL query.


texts: Trilateration using known points, but relative distances

texts:
I want to find the latitude and longitude of a point, based on knowing three nearby locations, similar to this well written post:

Trilateration using 3 latitude and longitude points, and 3 distances

The quirk in this situation is that I only know relative distances.  So for distances to three points (d1,d2,d3) I would only know that:

d1 = 4 * d3
d2 = 2 * d3

Is there a formula that lets me find x1,y1 while minimizing d3?  As in, what is the coordinates of the trilateration point where d3 is the smallest value that solves for the equation above?


texts:
I'm using a free acount an QGIS Cloud and I've uploaded the files such as in this tutorial: https://qgiscloud.com/en/pages/quickstart.

The problem is that the legend shows the layers but the map is empty. What seems to be the problem?

The website is https://qgiscloud.com/cove_adrian/WebGIS_uat_Victoria/


texts: Compability problems with OpenGeo Suite and already installed GeoServer and others

texts:
Sorry I couldn´t find the answer anywhere else.

I already have GeoServer 2.3.2, Postgres 9.2, OpenLayers 2.13 installed but in addition I need GeoWebCache and probably GeoExt.

Is it wice to download now the whole OpenGeo Suite even if the versions are different or should I try to download GeoWebCache and GeoExt separately?

Thanks for help!

Laura


texts: Control points in IAdjustXForm for georeferencing of raster (ArcObjects)

texts:
I am trying to create a georeferencing tool using ArcObjects in VB.Net. Though only a little text is available on this topic in forums as well as from ESRI; I found an esri article that states the use of IAdjustXForm for performing the desired transformations on input raster map and this article to create the AdjustXForm using control points.

I have built the interface to place a selected raster file as IRasterLayer on Map Control without any projection. After loading, user can mark the corners of the map image (un-referenced) on map control (by mouse click) and input the correct lat/lon for those points.

Dim sourcePC As IPointCollection = New Multipoint
sourcePC.AddPoint(TopLeftPoint)
sourcePC.AddPoint(BottomRightPoint)

I have tried building two IPointCollection initialized as a MultiPoint object by adding the user selected points (from mouse clicks) as source points (IPoint) and input coordinates (from text blocks) as target points (IPoint).

But it throws an exception - Exception from HRESULT: 0x80041034

So, now question is - I am not sure that control points for IAdjustXForm are meant to be the ones I have selected, or if yes, Are the source points pixel locations in image space or geographic coordinates after placing the raster on map?


texts: How can I draw features on GeoTIFF in Python?

texts:
I have a set of data (coords + value) and want to turn it into a heatmap in GeoTIFF. Should I calculate pixels and set them as array, or there are shortcuts?

This is what I do now:

    driver = gdal.GetDriverByName('GTiff')
    dst_ds = driver.Create('%s-%s.tif' % (name, zoom), height, width, 1, gdal.GDT_Float32)
    raster = numpy.zeros((width, height))

    cur.execute('select * from image')
    for i, row in enumerate(cur):
        raster[row[0], row[1]] = row[2]

I'd like to be able to do something like this:

 band = dst_ds.GetRasterBand(1)
 for row in cursor:
     band.drawCircle((row[0], row[1]), 20, row[2])

The important part is saving the depth of data.


texts: Dynamic selection of layer in QGIS atlas

texts:
I am making an atlas where I want to expose features with specific values. A good way to illustrate this could be this real estate database:

01 estate_a expensive yellowhouse
02 estate_b expensive yellowhouse
03 estate_c cheap     greenhouse
04 estate_d cheap     brownhouse
05 estate_e expensive greenhouse

The current atlas functionality in QGIS (2.18.0) supports making an atlas with one map for each estate. I want to make one map with all expensive estates and one with all cheap estates. Is this possible?

Also - can I have different symbols according to different house colors?


texts: Vector averaging wind direction using focal buffers in R?

texts:
I have a raster of 1 km wind direction data that I need to run a focal buffer on. Instead of using the default buffering mean function, I need to use vector averaging since my wind direction is in degrees (e.g. (360+0)/2 = 180/South, instead of 360/North!).

I tried using the circular mean function from the circular package, but this doesn't work with rasters.

How do I perform vector average focal buffers on rasters?

I wrapped avgWind function around circular mean but still get the same error.

library(raster)
library(circular)

# creating a sample raster for stackoverflow
xy <- matrix(sample(0:359), 10, 10)
r <- raster(xy)

avgWind <- function(r, units) {
  r_avg <- mean(circular(r, units = "degrees"))
  return(r_avg)
}

r_buf <- focal(r, focalWeight(r, d=20, type='circle')
               ,fun = avgWind(r)
               ,na.rm = T)

Error in x/180 : non-numeric argument to binary operator
In addition: Warning message:
In class(x) <- c("circular", cl) :
  Setting class(x) to multiple strings ("circular", "RasterLayer", ...); result will no longer be an S4 object

EDIT

After applying obrl_soil's example to my projected data the buffered data appears to be inverted.


texts: I have found that to carry out a kernel density estimate others are using the plugin AniMove. Please could someone point me in the right direction as I have looked in the Processing options and AniMove is not there under the providers section.


texts: How to save a google map direction into a postgresql/postgis database?

texts:
I have taken a direction(route) from google map.

Can I save this resultant direction/path/route to my postgresql/postgis database? How?


texts: Can a folder connection in ArcMap Catalog window be added using Python?

texts: I have created a set of general scripting tools for map automation, etc. In one of my tools a geodatabase is created in the workspace folder the script points to. Just out of curiosity, is there any way I can add a folder connection in the ArcMap catalog window once the GDB is created so the user doesn't have to manually create one (assuming they haven't already done so), or is that out of ArcPy's scope?


texts: What does Zero pixel value in NDBI mean?

texts: NDBI > 0 indicates built up, while NDVI <0 indicates no built up value. What does NDBI = 0 mean?


texts: I'm missing the Editor Option for Code Style SQL. Its not in the Options and I am using PyCharm Community 2019.3 on Manjaro.


texts: I have an issue where leaflet when embedded in Drupal will not allow the user to click on the map as the Y coordinates have a decimal on for some reason. As a quick fix, I am thinking to put the map into full screen mode outside of Drupal. However I would like to place a link/button on the map that points back to my Drupal site. Any suggestions on placing custom buttons on the map?


texts: From time to time it happens to me that rasters disappear at certain zoom level in ArcMap, or only some parts of the raster (strips) show up. I'm working on WinXP, ArcGIS / ArcView 10, Sp2, rasters are jpg files, both with and without pyramids. Anyone facing a similar behaviour? Any tips what is the reason and solution for this? It's really annoying.


texts: Cannot get correct Lat & Long scales on outside of OS UK map

texts: Using OS data and CRS of EPSG : 27700 on QGIS, on the new Print Composer how do I get the correct Longitude and latitude scale on the outside of the map? at the moment both scales on the X:Y axis start at 0, not as the actual map states?


texts: GRASS r.grow doesn't seem to change raster values at all

texts:
I have a raster map in GRASS that looks like this:

The red cells have a value of 1, the yellow cells have a value of 0.
I am trying to use the r.grow command to 'fill in' the gaps in the red lines, but it seems that whatever options I give the command, the output is either exactly the same as the input, or entirely 1, or entirely 0.
The way to run the command that seems most sensible is:
r.grow input=extracted@PERMANENT output=gr8 radius=3.01

which, as far as I can tell, should grow all positive cells by 3 cells. However, it doesn't seem to do anything at all - the output map is exactly the same as the input map. Changing the old and new values doesn't seem to change it either.
Am I doing something wrong here? I have read through the manual page for the r.grow command many times, but can't see what I'm missing.


texts: Adding attribute list to combobox in QGIS plugin?

texts:
I am developing a simple plugin at the moment which contains 3 combo boxes:

comboBox_1
comboBox_2
comboBox_3

Combobox 1 & 2 lists the vector layers available, however I would like the 3rd combobox to list the attribute column that both layers have in common (from the attribute table).

    # run method that performs all the real work
def run(self):
    # show the dialog
    self.dlg.show()
    for layer in self.iface.legendInterface().layers():
        if layer.type() == QgsMapLayer.VectorLayer:
            self.dlg.comboBox.addItem(layer.name())
            self.dlg.comboBox_2.addItem(layer.name())


texts: Linestring moving to the other end of the map near the edges


texts: I have a large polygon with many vertices that has a hole in its center. I want a polygon of the hole, the empty space. How can I do that, as it is much too complicated to digitize manually?


texts: I have GPS data of a route which is about 32 km. The number of recorded gps waypoints for this route is 1670. But I would like to have a waypoint every 10 meters. So the number of waypoints for this route should be increased to 3200.


texts: Importing Leica total station .GSI file data into ArcGIS Desktop 9.x?

texts:
I have a total station data file (test.gsi) that was collected with a Leica 1200.

How can I import the data to ArcGIS 9.x and plot points with the point name given by the total station observer. After importing the points i need to draw polygons around the points to define land parcels.


texts: Combining touching polygons, keeping all attributes from Class "A" feature, using ArcGIS Desktop?

texts:
How to automatically merge polygons sharing a common boundary, keeping all attributes from class "A" features?

Interactively I would do this in an edit session by:

selecting all objects of Name "One",
Editor toolbar >> Merge,
choose Class "A" to be the parent/master feature.

Repeat for features of Name "Two", "Three", ...

Dissolve solution from How to combine adjacent polygons sharing similar trait into single polygon? doesn't work because we need to keep other attributes which differ. In other words:

Dissolve/merge touching polygons where geometry [Name] is identical
Keep attributes from polygon with [Class] A (if multiple Class A's, ok to just keep first one)
Discard attributes from other features


texts: Ordnance Survey grids - download or generate?

texts: Does anyone know a source from which to download vector UK Ordnance Survey grids at varying resolutions? Or failing that, a tool that can create an OS grid at a specified resolution? There are plenty of ways to generate a vector grid, but I've found no way of populating the attributes with OS grid references.


texts: Converting geometry_dump to Point geometry vertex

texts:
If this was a Point geometry 010100000001DE02098A1B52C090A0F831E62E4640, is it possible to
represent the above geometry_dump Point Vertex. I am trying with St_collection...,


texts:
I'm trying to work out least cost paths in QGIS. I have a DEM and a slope raster extracted from it, and a set of points, between which I want to find the least cost paths. I can find a SAGA plugin, but it needs a cost raster.

How do I go about making one of those, and then getting the least cost paths?


texts: Masking Labeles in the Compositor in QGIS 2.8.1

texts:
I am looking for a way to remove/mask out Labels that are under a composer Element (Overview Map), which I moved in front of my main map.
Is there no way to do this? Maybe using the ItemID?
I tried the mask add-on, but it only seems to work with Polygons, which does not help me here, because the maps are generated by the Atlas function.


texts: QGIS 2.14 How to save and load parameter values in Batch Processing Tool

texts:
Am relatively new to QGIS, I am trying to save a query I carry out frequently using the Batch Processing - Extract by location tool. For example with the following parameters:

I then want to save these parameters so I can use them next time rather than add in and tick all the boxes etc etc because all the parameters stay the same pretty much. So I go to the Save function under parameters and save. But when I go to reload  the layers all appear with quotation marks which makes running the search invalid. What am I doing wrong?


texts: Annotating images on QGIS or using widgets?

texts:
I've georeferenced an old map on qgis and layered it on top of OSM, I want to do a class project and for it I would like to pinpoint certain areas on the map of which I have old photographs. I would like to know if it's possible to pinpoint the locations and attach photographs, so you can see photos of what that place used to look like, as well as write textual information. I know that you can do text annotations, but how about images?

I would like it to be something similar to like when you get locations pinpointed on google maps and when you press on the pinpoint you get the photo and information.


texts: Using points of elevation stored as Google Fusion Table to draw elevation profile

texts:
I recently extracted a huge number of elevation data from government plans and has inputted to Google Fusion Table. Now I wish to obtain a elevation profile from these data along a route, e.g. a bus route. Is there any method, e.g. using GIS software, to accomplish the above mission?

Please see my sample of elevation data collected in the format of Google Fusion Table here. If it is advised to convert the data into other format I will do so.
https://www.google.com/fusiontables/DataSource?docid=1tBcHh6y1Ri62zAbZnFo9jr0NQlLZkLK6VuD5Btd5

BTW I am just a newbie to GIS.


texts: Esri/ArcGIS Database Connection IWorkspaceFactory.Open COMException HRESULT 0x8004021A -2147220966

texts:
I am having a problem in connecting to the database through code in VS 2010 (C#). My colleagues don't have this problem on their machines and can login without any problem.

ArcGIS version 10.2.2, the same as my colleagues.
I can connect to the database outside the ArcMap, without any problem.

Also, I didn't find anyone with the same error code as mine "HRESULT 0x8004021A or -2147220966" in the Internet and already tried to search Google page 2 and forward so I'm kind of desperate and running out of ideas.

Here are some prints of the issue:

The IWorkspaceFactory.Open pops-up the ArcGIS Database Connection, and after I try to connect (which seems to not even trying to connect since its too fast in showing the error message) it just shows that error as you can see in the print. I can't see the IWorkspaceFactory.Open method since its a black boxed API. I already reinstalled the ArcObjectsSDKNet because the problem seems to be in the API.


texts: Is spat the correct flag to use in my ogr2ogr command?

texts:
I have an ogr2ogr command that I'm using to convert part of a Postgres database to a Geojson. I ingested all of Canada in my database and to start, I want to grab all of the buildings in the Toronto area for my first Geojson.

I've tried using spat in my command, but I don't get any features back. Is spat the correct flag or do I need to be using another one?

Command:

ogr2ogr -f "GeoJSON" ./toronto-buildings.geojson -t_srs 'EPSG:4326' PG:"host=localhost port=5432 user='MY USER NAME' dbname='gis'" -spat -79.9352, 43.4968, -79.1414, 43.9124 -sql 'SELECT way, "building" AS type FROM planet_osm_polygon ORDER BY ST_YMin(ST_Envelope(way)) DESC' -progress

Output:

{
   "type":"FeatureCollection",
   "crs":{
      "type":"name",
      "properties":{
         "name":"urn:ogc:def:crs:OGC:1.3:CRS84"
      }
   },
   "features":[

   ]
}


texts: Find nearby cities by a driving distance radius

texts:
I'm studying GIS services to solve one specific need:

Given x/y coordinates of one city, I want a JSON result that lists only cities (not facilities) that are nearby the specified city, within a driving distance radius (e.g. 100km).

I came close to that with two: Geonames and ArcGIS, but I couldn't extract this kind of data on neither (maybe it's my lack of experience).

Do you know a service that I can use for this specific need? (I don't want map rendering). If possible, can you give a query example using it?


texts: label placment in multipolygon with cartocss (Denmark, Greenland)

texts:
I am rendering an interactive world map with Tilemill and I have the problem that the country name for Denmark is placed onto Greenland. This is because Denmark is represented as a multipolygon with polygons for Denmark proper and Greenland, as Greenland is a territorial dependency of Denmark, and because Tilemill only places one text label onto a multipolygon .

My question is if it is possible to influence with Tilemill's cartocss on which of a multipolygon's subpolygons the textlabel is placed. I'd prefer for the text label to appear on continental Denmark.

I have tried the

    text-placement-type: simple;
    text-placements: 'E,W';

properties, hoping that this would select the more east or westward polygon, but to no avail.
I have also looked for a way to position a textlabel at an absolute position, but this also does not seem to be possible.

Finally, I would like to avoid manipulating the shapefile, as I am taking it from an external source.

Thanks!


texts: Gaussian plume model in PostGIS or geocouch

texts:
Does anyone know of existing code for dispersion modeling using Gaussian Plume or comparable model, that I could plug in / adapt to be used programmatically from PostGIS or alternatively via GeoJSON/geocouch?

I'm trying to model something like this: http://ready.arl.noaa.gov/gaussian-bin/dispsrc.pl to calculate dispersion values providing as few inputs as reasonably possible: lat/lon, Wind Speed, Wind Direction (possibly cloud cover  maybe elevation, cloud cover, etc)

My primary objective is to return computed values for % concentrations at points A, B, C, D downwind from emission point X, assuming 100% at point of release X

Secondary objective is visualize it along the lines of something like below


texts: Should I delete the existing replica when create a new replica with the same name?

texts:
I have successfully create a replica between 2 geodatabase and use a scheduler to synchronize the replica between the 2 geodatabase several months ago. But the scheduler didn't work because of changes in the original geodb. In this case I need to create the replica once again (according to user's request), my question is if I use the same replica name, is there a need of deleting the existing one? What will happen if I create the replica directly? Any potential dangers? Since this is a in operation system, I cannot afford of failure.

EDIT: After my research from internet, I even cannot find a way to delete a replica.


texts: Retaining last position of vessel position in geodataframe including time series?

texts:
I have geodataframe about vessel position which include unique code called 'mmsi' and time series. In the following 'mmsi' column has more than 1 'mmsi' code and time series (* mark).

           mmsi                 time
98    525107006  2018-12-07T23:55:48
295   563038200  2018-12-07T23:52:42
627   477162400  2018-12-07T23:49:50
704   244850097  2018-12-07T23:47:41
1141  352366000  2018-12-07T23:30:00*
1166  376265000  2018-12-07T23:28:29
1319  352366000  2018-12-07T23:19:53*
1345  477858800  2018-12-07T23:18:14

How to retain last position of a vessel of 1 'mmsi' considering the latest time and remove the rest?

I expect that it can be execute in a loop for other 'mmsi' code that have same issue.


texts: Orfeo Toolbox output rasters are 300x larger than inputs of the same size

texts:
I am trying to use Orfeo Toolbox CLI to do some processing, e.g. Smoothing:
otbcli_Smoothing -in path/to/mosaic.tif -out output/smoothing.tif uint8 -type mean

For now, the input, mosaic.tif, is just full of zeros and datatype uint8. The problem is that while the input file is approx. 9Kb, the output is 3.7Mb.
When I run gdalinfo mosaic.tif:
> Driver: GTiff/GeoTIFF Files:
> Size is 1014, 722 Coordinate System is: PROJCRS["ETRS89 / TM35FIN(E,N)",
> ...
> Image Structure Metadata:  COMPRESSION=LZW  INTERLEAVE=PIXEL
> ...
>
> Band 1 Block=256x256 Type=Byte, ColorInterp=Gray
> Band 2 Block=256x256 Type=Byte, ColorInterp=Undefined
> Band 3 Block=256x256 Type=Byte, ColorInterp=Undefined
> Band 4 Block=256x256 Type=Byte, ColorInterp=Undefined
> Band 5 Block=256x256 Type=Byte, ColorInterp=Undefined

When I run gdalinfo smoothing.tif:
> Driver: GTiff/GeoTIFF Files:
> Size is 1014, 722 Coordinate System is: PROJCRS["ETRS89 / TM35FIN(E,N)",
> ...
> Image Structure Metadata:   INTERLEAVE=PIXEL
> ...
> Band 1 Block=1014x1 Type=Byte, ColorInterp=Gray
> Band 2 Block=1014x1 Type=Byte, ColorInterp=Undefined
> Band 3 Block=1014x1 Type=Byte, ColorInterp=Undefined
> Band 4 Block=1014x1 Type=Byte, ColorInterp=Undefined
> Band 5 Block=1014x1 Type=Byte, ColorInterp=Undefined

So there are some differences (the part not shown is identical) but I don't know if (and how) that explains the explosion of the file size. I've checked that the output is the same datatype, uint8.
Question: What causes the output to be so much larger and how can I deal with it?


texts:
I have two shapefiles and I would like to know the distance between each point on File1 and each corresponding point on File2 (one to one relationship).

I have been using the Ruler and manually calculating this, but is there a more efficient way?

I am able to produce one table with the lat/long from Table 1 and the lat/long from Table 2. However, I do not know how to populate a distance column with an equation (and I do not know Python).


texts:
I'm downscaling some climate data for a small area in Italy, and have created the 36 monthly files for Tmin, Tmax and Prec. Having exported from ArcGIS 10.4 into ascii format (integers), then imported in DIVA-GIS 7.5 with no problems, the creation of the CLM file is hitting a snag. I'm getting the error message "attempting to read beyond end of grid file", then the module fails with an access violation.

I've tried running the program as Administrator to get around any potential read/write issues. I've previously successfully created CLM files with 30 arc second WorldClim data. Has anyone encountered this before?


texts: Geostatistical layer extrapolation in ArcGIS Pro

texts:
Just recently switched to using ArcGIS Pro from ArcMap 10.6
Does anyone know how to change the extent of a geostatistical layer in Pro?

In ArcMap this is a simple task (as outlined here) however this option does not appear in the layer properties in Pro.


texts: projection/pyproj puzzle, and understanding SRS format

texts:
I have a shapefile that has a .prj file, with the following wkt contents:

PROJCS["NAD_1927_StatePlane_New_York_West_FIPS_3103",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-78.58333333333333],PARAMETER["Scale_Factor",0.9999375],PARAMETER["Latitude_Of_Origin",40.0],UNIT["Foot_US",0.3048006096012192]]

I am using the following python code to convert that to a pyproj projection:

prjText = open( prjPath, 'r').read()
srs = osr.SpatialReference()
if ( srs.ImportFromWkt( prjText ) ):
    return ( False, "error importing .prj information from " + prjPath )
inProjection = pyproj.Proj( srs.ExportToProj4() )

I then create another projection:

outProjection = pyproj.Proj( "+proj=latlong +datum=WGS84" )

and then go through each point p from the shapefile and convert it:

print "before transform: " + str( p[ 0 ] ) + ", " + str( p[ 1 ] )
p = pyproj.transform( inProjection, outProjection, p[ 0 ], p[ 1 ] )
print "    -> after transform: " + str( p[ 0 ] ) + ", " + str( p[ 1 ] )

when I print the srs.ExportToProj4() text, it gives:

+proj=tmerc +lat_0=40 +lon_0=-78.58333333333333 +k=0.9999375 +x_0=152400.3048006096 +y_0=0 +ellps=clrk66 +to_meter=0.3048006096012192 +no_defs

Furthermore, using another program that I have, I import the same shapefile without using the .prj file. To do this, I specify that the projection is "State Plane", and specifically that it is "New York West NAD27 (ftUS)". This other program then determines that that means the same as epsg:32017, and indeed when I look up epsg:32017 I see that its srs string matches the one above.

Thus, it seems like my inProjection variable (above) should be set up correctly. However, when I run the code above to transform each of the points, I get, for the first point:

before transform: 759232.003438, 1149854.52147
    -> after transform: -70.1071453809, 50.036825617

The puzzle:

That lon/lat, approximately (-70, 50), is not in the right place. By importing the file into either the above-mentioned program, or into ArcGIS Explorer Desktop, I see the point (correctly) in northern New York, just south of Lake Ontario. But (-70, 50) is up in central Quebec. So it's off by like 1000km to the north and east.

So the questions are:

(1) Does it seems likely something is going wrong in pyproj (I'm using version 1.8.9)? (since according to these other programs my srs string is correct, so everything is good up to the point I use pyproj)

(2) Based on the wkt and srs string shown above, what lon/lat should I expect the the shapefile point (759232.003438, 1149854.52147) to correspond to? Is that point the x/y in meters relative to the lon_0 and lat_0 of the srs string? If that were the case, I would indeed expect a point somewhere up in Quebec, and not in northern New York. However, since both of these other programs are correctly mapping the point to northern New York, that must not be the correct interpretation of the file coordinates with respect to the file's projection?

Thanks for any help.


texts: Convert from GDA2020 to GDA94 using proj string in PostGIS

texts:
Following this answer, I was able to convert from GDA94 to GDA2020 using proj string (I'm using PostGIS).

Now is there a way to convert it back? I tried to invert the signs from @JasonD answer but it didn't work:

+proj=longlat +ellps=GRS80 +towgs84=0.06155,-0.01087,-0.04019,0.0394924,0.0327221,0.0328979

When I use the above string on the converted geometry the coordinates are not the same from what they were before the conversion.

Can someone advise how to do it?


texts: Error when creating spatial index in MongoDB

texts:
I am wanting to use MongoDB for some basic spatial calculations (on some fairly large data), but I get an error when trying to create an index on the US Counties shapefile retrieved from the US Census' TIGER data. I had to do some manipulation in the terminal to get this data into MongoDB in the first place since Mongo can't ingest shapefiles directly:

$ ogr2ogr -f GeoJSON counties.geojson tl_2015_us_county.shp -t_srs EPSG:4326 # covert shapefile to geojson
$ jq --compact-output ".features" counties.geojson > counties_reformatted.geojson # reformat the .geojson because mongoimport throws an error if you don't
$ mongoimport -d mydb -c counties < counties_reformatted.geojson --jsonArray --batchSize 1 --drop # import to mongodb

I then create the index in the mongo shell with

> db.counties.createIndex( { "geometry" : "2dsphere" } )

but I get this error:

.
.
.
Edges 839 and 841 cross. Edge locations in degrees: [-99.8924200, 36.5932380]-[-99.8960610, 36.5932360] and [-99.8960600, 36.5932360]-[-99.8960630, 36.5932360]",
    "code" : 16755

This is along the border of Oklahoma and Kansas.

I have been able to get this process to work with shapefiles of Mexico's states, so it has worked in the past. Where is the error coming from? The ogr2ogr conversion? Or something wrong with the shapefile in the first place?

The counties, tl_2015_us_county.shp have EPSG:4269 by default, so I suspected there may be a problem in converting to EPSG:4326, but changing t_srs to EPSG:4269 in my ogr2ogr sciprt did not fix the problem. I need (I think) the counties in EPSG:4326 because I have a much, much bigger points layer already in this system.

Notes:

Maybe Stack Overflow would be a better place for this?

Inverted lat/lon is not the problem, and this workaround does not work for me.

Update:

I also tried using two different shapefiles (with different resolutions) from the US Census here, and I get a very similar error. When plotting the points that mongo says overlap, there is not a perfect coincidence with the points of the shapefile. This leads me to think that something is going wrong in my conversion process.

I'm guessing this could be solved by using really low resolution data, but I'm not really comfortable with that.

Also, spatial queries without the index are not really an option here. They take excruciatingly long without the index.

Version info:

os: Ubuntu 14.04
mongodb: 3.2.7
ogr2ogr4: GDAL 1.10.1, released 2013/08/26
jq: 1.3

Shapefile (in case anyone wants to replicate this):

wget ftp://ftp2.census.gov/geo/tiger/TIGER2015/COUNTY/tl_2015_us_county.zip --no-parent --relative --recursive --level=2 --accept=zip --mirror